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FUNDAMENTALS  OF  THE  EPCE-FUNCTION  METHOD  VIA  LAPLACE'S  EQUATION 


by 

P.  M.  QUINLAN^ 


Introduction 

1.  OUTLINE 

The  Edge-Function  Method  as  applied  to  two-dimensional  problems  is  a 
procedure  for  obtaining  approximate  solutions  to  boundary  value  problems  in 
regular  or  irregular  polygonal  regions  with,  or  without,  cavities  and  cracks. 
It  may  be  described  as  a piecing  together  of  "asymptotic"  solutions  to  a set 
of  linear  partial  differential  equations  for  the  several  parts  of  a domain 
D to  satisfy  the  boundary  conditions  in  a discrete  least  squares  sense. 

The  method  was  originally  developed  by  Quinlan  [1]  for  the  solution  of 
the  torsion  problem  for  prismatic  bars  of  polygonal  cross-secticm.  In  the 
succeeding  10  years  it  has  been  successfully  applied  to  problems  in  the  bend- 
ing of  isotropic  thin  plates  [2,3,4];  coupled  linear  systems  in  elastostatic 
[4, 5, 6, 7]  and  in  moderately  thick  plates  and  shallow  shells  [8];  cracks  and 
stress  concentrations  in  elastostatics  [9];  vibrations  of  Thin  Plates  [10] 
and  vibrations  of  Shallow  Shells  [11], 

The  present  paper  arose  out  of  a course  of  lectures  on  E.F.M.  given  by 
the  author  in  December  1975  at  The  Georgia  Institute  of  Technology.  The  aim 
is  to  illustrate  as  simply  as  possible  the  main  algebraic  and  programming 
features  of  the  method  by  applying  it  to  Laplace's  equation  as  it  arises  in 
heat  flow  and  torsion  problems. 

2.  PLAN  OF  PAPER 

The  paper  is  divided  into  four  sections: 

(a)  Section  1:  Introduces  Edge-Functions  and  Vertex  Functions  for  steady 

state  heat  flow  in  a polygonal  region  the  temperatures  being  specified 
on  the  boundary.  The  resulting  boundary  identities  are  satisfied  by 
harmonic  matching  using  Fourier  sine  series. 

The  importance  of  the  vertex  equations  is  shown.  The  Edge- 
Functions  are  obtained  using  a half-strip  mathematical  model. 

A simple  program,  LAPEX,  and  an  illustrative  example,  is  given  in 
Appendix  A for  problems  covered  by  this  section. 


(1)  Professor  of  Mathematical  Physics,  University  College,  Cork,  Ireland. 


(b)  Section  2:  Extension  is  made  to  regions  with  mixed  conditions  on 

the  boundary  segments, and  circular  holes  are  included.  Harmonic 
matching  is  now  based  on  using  the  (full)  Fourier  series,  with 
corresponding  (full)  Edge-Functions. 

Provision  is  made  for  Logarithmic  Singularities,  where  they 
arise,  using  Log-Vertex  Functions. 

Programming  requirements  for  a comprehensive  program  LAPGEN, 
as  a development  of  LAPEX,  are  discussed  and  a suitable  program 
LAPGEN,  with  examples,  is  attached  in  Appendix  B. 

(c)  Section  3 - Curved  Boundaries:  This  uses  conformal  mapping  to  extend 

the  treatment  of  the  previous  section  to  holes  of  elliptic  shape,  and 
to  circular  and/or  elliptic  indentations  on  the  boundary.  The  ensuing 
Curved  Edge-Functions  are  developed,  and  programmed  in  subroutine  PMAP 
in  LAPGEN  for  elliptic  curves. 

(d)  Section  4 - Other  Basic  Problems:  This  section  treats  the  case  of 

a reentrant  angle,  and  also  where  singularities  occur  on  the  boundary, 
using  appropriate  half-strip  models.  These  features  have  not  been 
included  in  LAPGEN,  but  are  available  in  QUINP. 

The  models  used  may  be  termed  "BASIC  MODELS"  for  the  problem, 
the  corresponding  regions  being  termed  "BASIC  REGIONS"  and  the 
corresponding  functions  being  termed  "BASIC  FUNCTIONS". 


(e)  Section  5 - Harmonic  Fitting:  A discrete  least  squares  method  fotf 

fitting  a set  of  orthogonal  functions  to  a given  curve  is  developed. 
It  is  applied  to  solve  the  boundary  identity  problem  by  minimising 
the  boundary  residuals,  on  each  segment,  using  a discrete  least 
squares  criterion.  This  process  is  given  the  name  "Harmonic 
Fitting" , and  follows  from  Harmonic  Matching  when  the  integrals 
concerned  are  replaced  by  the  corresponding  trapezoidal  rule 
quadrature  formula. 

(f)  Section  6 - Distinctive  Features  of  the  Edge  Function  Method. 

Appendix  A:  ' LAPEX,  with  example  of  heat  flow  in  a quadrilateral. 


Appendix  B: 
Appendix  C* 


LAPGEN  with  example  involving  the  torsion  of  a quadilateiral 
region  with  an  elliptic  hole 

Data  Input  for  LAPGEN  on  the  "black  box"  principle. 
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SUMMARY  OF  METHOD 

Accordingly  to  apply  The  Edge-Function  Method  to  any  new  situation 
e.g.  thin  plates,  shells,  3-D  Elasticity,  we  require  to  find  the  corresponding 
BASIC  FUNCTIONS,  and  fit  a superposition  of  these  to  the  boundary  conditions 
by  "HARMONIC  FITTING"  using  a computation  scheme  similar  to  that  given  in 
the  attached  program  LAPGEN . 

Sections  3,  4 and  5 may  be  omitted  when  studying  the  Edge-Function  Method 
for  the  first  time. 


SECTION  1 


> x 


Consider  the  linear  partial  differential  equation: 

,2 

L(u)  = f(x,y); 


L = +jl_ 

3x2  3y2 


(1) 


where  L,  in  the  present  paper,  is  restricted  to  the  Laplace  operator. 

It  is  required  to  find  a function  u = u(x,y)  which  satisfies  the 
above  equation  in  a domainJ^__bounded  by  any  number  q'  piecewise  smooth 
segments,  or  straight  lines,  and  which  takes  prescribed  values 
u = gq(x,y) 

on  the  various  boundary  segments  corresponding  to  q = l,2,...,q’. 

For  example  in  a steady  state  heat  flow  problem  u is  the 
temperature  and  f(x,y)  is  the  heat  generation  term,  and  for  a heat 
source  of  strength  Q at  a point  (u,v),  f is  given  by 
f(x,y)  = Q5(x-u,  y-v), 
where  6 is  the  delta  function. 

If  u = u**  is  a particular  solution  of  eqn.(l),  then 
L(uP)  = f (x,y) 

A more  general  solution,  involving  the  addition  of  another 

function  u , may  be  written  in  the  form 

p . c 

u = ur  + u 

where,  since  L is  linear, 

L(uP  + uC)  = L(uP)  + L(uC)  = f(x,y), 

and  on  subtracting  eqn.(4),  it  follows  that 

L(uC)  = 0, 

which  is  called  the  complementary  equation  of  eqn.(l). 


(2) 


(3) 

(4) 

(5) 

(6) 
(7) 


The  problem  is  now  split  into  finding 

(a)  a particular  solution  u*5,  and 

(b)  suitable  complementary  functions  u , with  sufficient  flexibility 
to  enable  solution  u1^  + u to  satisfy,  in  an  approximate  manner, 
the  boundary  conditions  (2) . 

The  present  paper  is  restricted  to  the  case  f = 0, 
for  which  u^  = 0. 

On  obtaining  I suitable  solutions  u? , i = 1,1;  of  eqn.(7) 
and  superposing,  we  obtain 


involving  the  I arbitrary  superposition  constants  A^,  and  these  must 
be  evaluated  to  satisfy  conditions  (.)  in  an  acceptable  approximate 


manner. 


I eqns.  are  required  to  determine  the  I constants  A^.  Various  schemes 
for  producing  I eqns.  are: 

) Point-Matching,  or  collocation  [Conway,  Leissa  etc.] 

Select  I points  on  boundary  in  same  manner  - equidistant 
distribution,  or  have  a somewhat  closer  spacing  near  corners— 
and  write  eqns.  to  ensure  that  boundary  conditions  (2)  are 
satisfied  at  the  selected  points.  Check  the  boundary 
residuals  (or  deviations)  at  points  in  between  the  selected 
points.  The  system  of  simultaneous  equations  may  become 
unstable,  giving  ill-conditioned  equations,  as  I is  increased. 

There  is  an  art  in  selecting  the  optimum  spacing  for  the 
points. 

) Least  squares  fitting  This  requires  that  the  boundary  conditions 
be  satisfied,  in  a least  squares  sense,  at  M points, 
where  M > I,  there  being  more  points  selected  than  there 
are  unknowns.  This  improves  the  process  up  to  about  M = 21, 
but  much  additional  computation  is  involved  in  producing  the 
resulting  normalised  equations. 

) Harmonic  Matching  (Quinlan  1962-70) 

Consider  the  boundary  identity  on  the  q side,  in  fig.l, 
where  points  on  the  side  are  specified  by  the  parameter  x^, 
which  rcquireSjOn  substituting  from  eqn.(8)  in  (2)^  that 


1 Aiui  H sn(x,y), 

i-l  A L 4 


<£ 


for  all  points  P(xq)  on  q 


th 


side  in  the  range  (0,aq), 


Since  the  coordinates  (x,y)  for  any  point  P can  be  expressed 
in  terms  of  xq,  identity  (9)  can  be  expressed  as 

I A.uJ(x')  - gq(x’)  = 0, 


i=l 


or,  for  short. 


*q<*^  H °i  0 aq» 


with  similar  identities  for  each  boundary  segment. 

On  expanding  identity  (11)  in  a Fourier  sine  series,  we 
obtain 


CO 

z 

N=1 


C„  sin  nx'  = 0; 
N q 


n = 


Nir 


where  the  Fourier  coefficients  are  given  by  the  integral 

, a 

C.,  =>  — f q ij)  (x')  sin  nx'  dx' 

N aq  o q q q q 

Since  we  have  but  I degrees  of  freedom  in  eqn.(8)  - or  I arbitrary 

constants  A.  - and  if  these  were  to  be  shared  cut  equally  between  • 

* 

the  q’  boundary  segments,  we  would  then  allocate  N.  , where 

N*  = I/q’ 

to  each  boundary  identity  (11).  Accordingly,  the  best  way 

* 

to  satisfy  identity  (12)  would  appear  to  be  to  set  the  first  N 
harmonics  to  zero,  giving  rise  to  the  equations 

ci  = 0 


(10) 


(11) 


(12) 


.(13) 


(14) 


(15) 


v° 

This  process  is  called  Harmonic  Matching 

On  substituting  from  eqs.(ll)  and  (10)  into  (13),  we  obtain, on  setting 


CN  = 0,  for  the  Nth  harmonic: 

la  a 

. Z A.  f q u?(x')  sin  nx'  dx’  - f q g (x')sin  nx'  dx'  = 0, 
i-1  1 o 1 q q q o q q q q 

where  N goes  from  1 to  N . On  evaluating,  a linear  equation  results 
in  the  unknowns  A^.  Similar- linear  equations  follow  for  each  side  of  the 
polygon  and  accordingly  a total  of  1 linear  equations  result  for  the  I 
unknowns,  A^,. which  system  may  be  written  in  the  matrix  form 

& - Q, 

where  P is  the  I x I coefficients  matrix. 


(16) 


(17) 


However  setting  the  first  N harmonics  to  zero  does  not 
satisfy  (12)  exactly,  but  leaves  as  a residual  the  higher  harmonics 


E C„,  sin  nx' : 
N=N  " q 


N = N +1 


and  we  require  this  to  be  negligible  - which  will  be  the  case 
only  if  the  series  is  sufficiently  convergent. 

As  shown  in  texts  on  Fourier  series  the  sine  series  expansion 
of  (11)  will  have  a Gibbs  effect  at  the  ends  x'  = 0;  x'  = a - 

1 q q q 

with  convergence  only  of  order  — - unless  the  function  being 

expanded,  <Jj_,  is  zero  at  both  ends  of  the  range.  In  this  latter 

q 3 

case  the  convergence  improves  considerably,  to  be  of  order  1/M  , 

and  consequently  the  residual  series  (18)  can  be  made  negligible, 

•jjf  * 

in  a practical  sense,  by  taking  N sufficiently  large.  Accordingly  we 
must  ensure  in  identity  (11)  that 


VXq>  = 0 ’ 


= 0 and  x'  = a 


This  involves  two  additional  equations  for  each  side  - called 

the  Vertex  Equations  - and  if  we  still  wish  to  restrict  to  N 

equations  per  side,  we  must  reduce  the  harmonic  equations  set 

(15)  for  each  side  by  two> corresponding  to  truncation  at  L harmonics 

where  L = N -2.  It  is  obvious,  due  to  the  greatly  increased 

convergence  of  the  expansion  (12)  to  1/M^,  that  this  scheme  is  a 

•k 

better  utilisation  of  the  available  resources  per  side  (N  constants) 
in  seeking  to  satisfy,  in  an  acceptable  approximate  way,  the 
boundary  identities  (11) . . ' 

If  the  temperature  is  continuous  at  all  points  on  the  boundary, 
including  all  corner  points,  then  eqs.(19)  reduce  to  a single  equation 
ij»  = 0 at  each  vertex. 


r 


Basic  Problems 

If  solution  form  (8)  is  capable  of  providing  an  acceptable 
solution  to  the  problem,  then  it  is  especially  important  that  it 
should  represent  the  main  features  of  the  solution  in  all  the 
critical  regions  of  Fig.  1.  On  examining  some  critical  parts  of 
the  boundary  we  note  that: 

(a)  in  any  corner,  say  at  A,  the  solution  (1-8)  must  approach 
the  solution  in  an  infinite  sector  B'AE',  where  B'  and  £' 
are  points  at  infinity  on  AB  and  AE  respectively.  In 
particular,  if  there  is  any  singular  behaviour  - physically 
corresponding  to  infinite  stresses  - in  the  infinite  sector 
solution,  there  must  be  corresponding  singularities  in 
solution  (1-8).  Hence  we  must  first  solve  the  infinite 
sector  problem,  and  include  its  more  singular  terms  - called 
Vertex  Functions  - in  superposition  (8)  for  each  of  the 
vertices  in  Figure  1. 

(b)  at  a distance  from  any  corner  where  the  vertex  effects  have 
moderated,  say  on  A*B*,  solution  (8)  should,  in  the  neighbour- 
hood of  A*B* , behave  like  the  solution  for  the  half-plane 
problem  bounded  by  A'B’,  where  A'  and  B'  are  points  at 
infinity  on  AB.  Such  solutions  are  termed  Edge-Functions . 
and  an  appropriate  number  of  these  must  be  included  for  each 
straight  sicle  of  Figure  1. 

(c)  in  the  vicinity  of  interior  hole  C the  solution  should 
approach  that  of  a hole  C in  an  infinite  region,  and  again  ' 
an  appropriate  number  of  the  more  singular  (or  characteristic) 
of  these  solutions  should  be  included. 

The  appropriate  "mix"  of  functions  from  (a) , (b)  and  (c)  will 
be  determined  later  from  analytic  considerations.  We  now  proceed  to 
solve  the  basic  problems  for  Laplace's  equation,  as  exemplified  by 
the  heat  flow  problem: - 


= 0, 


where  the  superscript  c has  been  omitted. 


t j yr.;[<n;x  FUNCTIONS  (Angle  Problem) 

Polar  type  solutions  of  Laplace's  equation  are  provided  by  the  well 
known  harmonics  of  order  X: 

u = ArX  cosXO  + BrX  sinX9, 

where  A and  B are  arbitrary  constants,  and  X is  arbitrary.  The 
appropriate  boundary  conditions  for  a sector  of  angle  a,  with  one 
side  9=0,  are  the  corresponding  zero  conditions, 
u = 0 ; 9 = 0 

u = 0 ; 9 = ctj 

as  these  give  rise  to  an  eigenvalue  problem. 

On  applying  conditions  (22)  to  eqn.(21a)  we  obtain 
0 = A 

0 = BrX  sinXa, 

a solution  to  which, other  than  the  trivial  one  u = 0,  is  given  by 
sinXa  = 0 

This  is  called  the  eigenvalue  equation  for  the  problem,  and  its 
solution  is 

Xa  = ksr, 

where  k is  any  integer.  On  denoting  the  above  discrete  set  of 
values  - called  eigenvalues  — for  X by  X^  we  obtain  as  solutions 

Xk 

u^  = Br  sinX^O;  X^  = kir/ct, 
where  u are  called  the  eigenfunctions  for  the  problem. 

K ^ 

Similar  eigenfunctions  - here  called  Vertex  Functions  - 
must  be  incorporated  in  "mix"  (8)  from  all  vertices  j = l,j';  ( j ’ =q ' ) • j 
Accordingly  the  kth  eigenfunction,  or  vertex  function,  for  the 
j*"*1  vertex  is  denoted  by  the  symbol  V,  . , where 


a, 

V.  . = u,  = B,  .[r  sinX,  9] . 
kj  k.  kj  k j 

the  subscript  j outside  the  square  bracket  denotes  that  the 

corresponding  origin  of  coordinates  is  at  vertex  j,  with  9=0 

along  side  j.  The  corresponding  arbitrary  constant  is  written 

as  B . , since  it  can  change  with  k and  jJ 
KJ 

Note  that  V.  . is  a point-function,  and  this  can  be  emphasised, 
kJ 

when  required,  by  writing  it  as  (P) . 

Eqn. (20)  is  invariant  under  a rotation  of  axes,  and  accordingly  on  taking 
the  polar  axis  9.=  0 along  the  jth  side  with  pole  at  vertex  j,  harmonic 
solutions  (21),  and  eigenfunctions  (24)  apply  to  the  vertex  j. 


10. 


The  permissable  range  of  X^ values  is  determined  by  the  physics 
of  the  problem.  Negative  values  of  X^  cannot  be  admitted  as  these 
would  require  infinite  temperatures  at  the  corresponding  vertices, 


but  can  infinite  temperature  derivatives  (at  r.  =0)  be  admitted? 
Consider  the  heat  flux  q,  where 


~ o(3u  A 1 dU  A, 

q - -BV»k  - -Btjjr  * vw9> 

The  flow  through  any  circular  arc  of  angle  a,  due  to  each 

3 V,  . 


V^j,is  given  by 


- ra 
"j  ' ;o 


3r 


rd0 


Xr  a 

= X^r  / sinX^0d0 
X 0 

= r {1  - cosX^a}, 

and  this  is  infinite  - and  thus  physically  inadmissable  - if 
X^  < 0.  Accordingly  in  (25)  we  require  that 

K 

The  corresponding  heat  flux,  as  given  by  (26)  on  substituting 

X. 


k-1 


, and  if  0 < X^  < 1 the  flux  is 


for  V,  . , is  of  order  r 
kj 

infinite  at  r =0,  but  this  is  admissable  since  the  nett  heat 
flow  through  the  vertex  is  zero. 

Note  that, if  X^  is  not  an  integer,  all  roots  X^  > 0 
have  infinite  derivatives  of  order  p and  higher, where  X^  - p < 0. 
Accordingly  all  functions  V,  . are  singularity  functions  in  the 
sense  that  all  their  derivatives  above  a certain  order  are 
singular  at  r = 0. 

When  a temperature  discontinuity  occurs  at  any  point  P on 

the  boundary,  or  at  any  vertex,  the  corresponding  discontinuity 

C3n  be  incorporated  in  "mix"  (8)  by  including  an  zero  order 

harmonic  solution  to  eqn.(20)  for  each  such  point  given  by 

u = B 0 
o o 


(26) 


(27) 


(28) 


(29) 


(b ] EDGE  FUNCTIONS  (half-plane  problem) 

Consider  the  solution  in  half  plane  y!  >0,  where  the  edge- 

th  * i 

axes  for  j side  denoted  (x! ,y!)  are  taken  respectively  along  and 


perpendicular  to  the  j 
reference  axis  (x,y). 


th 


J J 

side  and  are  at  an  angle  0^  to  the 


On  transforming  eqn.(7)  to  the  above  edge-axes  with  the  aid  of 
the  relations: 

x!  = (x-x.)cos(f,.  + (y-y.)sin^,. 


we  obtain 


where  the  coefficients  in  L'  may  depend  on  <f> _ . 

Since  Laplace's  equation  is  invariant  under  a rotation  of 
axes,  then  eqn. (7)  is  invariant  and  hence  transforms  to 


will  be  of  the 


Solutions  of  eqn. (32)  that  are  trigonometric  m x 
form 

u = F(y')sin(mx'  + a), 
and  on  substituting  in  (32)  it  follows  that 


the  solution  to  which  is 

F(yj)  = Ae"”myj  + Bemyj 

On  seeking  functions  that  could  represent  the  propogation  inwards 
of  the  decaying  effects  of  boundary  adjustments  on  side  j,  we  must 
exclude  the  positive  exponential  part.  Accordingly 


where  the  arbitrary  constant  A may  change  with  m and  j and  is 

written  as  A ••  We  define,  as  in  Quinlan  ( 1 ) , on  setting  u-0, 
mJ 

Edee-Functions  for  straight  edges  in  a Laplace  problem: 


where  further  analysis  is  required  to  show  that  m must  be  restricted 
to  the  discrete  set  (36). 


in  the  case  of  this  is  a point  function  and  can,  if  required, 

be  written  as  Kmj <p) • 

EDGE-FUNCTION'S  (half-strip  problem) 

Recent  developments  in  Edge-Functions  show  that  it  is  useful 

to  consider  the  basic  problem  (b)  as  approaching  the  solution  in  a 

semi-finite  strip  on  the  side  AB,  (y^  = 0)  as  base,  with  infinite 

sides  x'.  = 0 and  x!  = a..  Since  the  other  sides  of  the  boundary 
J J _ J _ * * 

are  assumed  to  have  little  effect  on  the  solution  along  A B , it 

follows  that  the  insertion  of  the  semi-infinite  sides,  with  appro- 
priate boundary  conditions,  will  likewise  have  little  effect  on 

& * 

the  solution  along  A B . 

Accordingly  in  the  case  of  sine  series  expansion  (12)  , we 
set  up  the  half-strip  problem,  for  Fig.  2 under,  with  the  boundary 
conditions 

u=0:x!=0,x!=a. 

J J J 

u ->  0 as  y\  ->  05 

J 

u = G(x!)  ; y!  = 0,  0 < x!  < a.. 

J j - J “ J 

We  seek  a solution,  in  the  usual  manner,  to  eqn.(32)  that  is 
trigonometric  in  x!  - corresponding  to  the  pair  of  zero  boundary 
conditions  - in  the  form 

u = F(yj)sin(tnxj  + a) 

where,  as  in  (33a),  on  solving  for  F(y!)  and  applying  condition 

u + 0 as  y!  + »,  we  obtain  ' 

J 

u = A . e sin(mx'.  + a) 
mj  J 

Conditions  u = 0 when  x!  = 0 and  x!  = a.  require 

J J J 

a = 0 ; m = irM/  a ^ ; H = 1 , 2 , . . . 

The  solution  for  an  actual  half-strip  follows,  in  the  usual  way,  on 
taking  the  superposition 

00  _ * 

. u = Z A . e ^ j sin  mx!  ; m = r M/a. 

M=1  mJ  3 3 

where  the  remaining  condition  (35)  requires  that 


G(x • ) = Z A . sin  mx!  ; 

3 M=1  ^ J 

the  coefficients  A . for  a Fourier  sine  series  expansion  being 

mj 

given  by  a. 

2 J 

A . = — | G(x!)sin  mx!  dx! 

“J  aj  o . J J J 


t 


V 


T„is  that  Edge- Functions  ofJaS-QZl. 

as  b"t  to  tlie  ■ 
"^^^77777  solution  W'  (81-  However  the  coeff  rcrcnt.sj^jsst 

77777777^0,00  scheme  V i.  that  aUous  for  the  efto^ 

7777  other  boundaries  - ralh.r  than  A„j  ..  Tven  >*  ««  <*»>* 

crude  mathematical  half-strip  model.  It  might  reasonably  be 

expected  that,  when  C(x')  is  not  identically  aero,  the  actual 

calculated  coefficients  would  be  of  the  same  order  as  these  result- 

ing  from  (38) . 

If  expansion  (12)  is  truncated  at  L harmonics,  it  follows 
that  the  series  (37)  should  be  truncated  similarly  at  M = L,  or 
in  solution  W (8)  , we  shouldj  nclude  as  many  Exjunctions 
from  each  side  as  we  have  si gni fi cant  - i.e.  retained_-_harmonics 


It  is  evident  that  there  is  a distinct  advantage  in  using  the  above 
half-strip  model  in  place  of  the  half-plane  model  in  20(b). 


It  remains  to  provide  additional  functions  equal  to  the 
total  number  of  vertex  eqs.(19)  in  the  problem  - and  the  obvious 
sources  to  provide  these  are  the  Vertex  Functions  V^. . It  has 
been  found  in  practice  that  it  does  not  matter  how  many  are 
taken  from  each  of  the  various  vertices,  provided  that  sufficient 
of  the  lower  eigenvalues  X^X^Xg...  are  included  from  each  vertex 
to  imbed  its  asymptotic  behaviour  into  the  solution  form. 

EXAMPLE 

It  is  a simple  programming  task  to  set  up  matrix  (17) > using 
eqs . (16)  and  (19)  with  truncation  at  any  specified  level  L.for  a 
polygonal  region  with  continuous  temperatures  on  the  boundary. 

A simple  program,  LAPEX,  is  given  in  Appendix  A(l) , with  some 
explanatory  captions. 

Appendix  A(2)  gives  results  for  a quadilateral  region  with  zero 
temperatures  on  three  edges,  while  on  the  remaining  edge  q = 1 

condition  (2)  is  specified  by 

u = g1(x,y)  = t(l  - t)  ; t = x[/al 

The  integrals  in  eqs . (16)  are  evaluated  by  the  trapezoidal 
rule,  noting  that  the  integrand  is  zero  at  both  ends  of  the  range 
since  sin  nx’  is  then  zero.  This  anticipates  the  replacement  of 
Harmonic  Matching  by  the  much  more  economical  Harmonic  Fitting 

Method  as  given  in  Section  5. 

The  equations  are  grouped  in  successive  harmonic  sets 

A _ j h = 2 N = L,  the  equations  from  each  side  being 

arranged  consecutively  in  each  set.  The  columns  for  the  Edge-  , 
Functions  and  the  rows  for  the  corresponding  harmonic  equations 
are  arranged  to  intersect  on  the  diagonal  of  the  coef f icieru-s 
matrix.  This  produces  diagonally  dominant  matrix  for  the  E.F. 
coefficients  - the  characteristic  matrix  form  for  E.F.M.  The  four 
vertex  equations  are  arranged  in  rows  after  the  harmonic  equations, 
and  the  columns  for  the  corresponding  vertex  functions  follow  those 
for  the  edge-functions.  The  right  hand  side  is  put  up  as  the  final 
column  of  the  matrix,  and  the  attached  simple  solution  subroutine, 
Qsolve,  solves  the  matrix  using  Gaussian  elimination. 

With  a view  to  their  extension  in  Appendix  B simple  subroutines 
EDGEF  and  POLW  are  introduced  for  the  evaluation  of  the  respective 
Edge  and  Vertex  Functions  at  any  point  (xk,yk). 


■zr*  j*. 


The  lie  si  duals  - the  difference  between  the  values 

computed  from  the  solution  form  (8)  as  obtained  and  the  prescribed 
boundary  values  - are  computed  and  printed  at  a specified  KDLV 
points  on  each  side.  Their  root  mean  square  is  obtained  and 
printed  as  a suitable  overall  measure  of  how  well  the  boundary 

conditions  are  satisfied  on  each  side. 

A single  data  card  specifying  points  Mu^v^  and  B(u2,v2> 
and  MDIV,  instructs  the  program  to  compute  the  temperature  at  each  o£  the 
MDIV  equidistant  points  on  the  line  AB.  In  the  production  phase 
use  is  made  of  that  part  of  the  program  set  up  for  the  calculatron 
of  the  boundary  residuals,  merely  by  cutting  the  j cycle  to 
single  pass,  corresponding  to  j - 1.  An  indicator  « code  - 1 
indicates  the  production  phase;  the  boundary  residuals  phase 
being  denoted  by  Ncode  - 0.  A blank  card,  which  reads  MDIV  as  zero. 

terminates  the  program. 


SECTION  2 

Mixed  Boundary  Conditions  and  Full  Harmonic  Matching 

We  now  proceed  to  generalise  to  mixed  boundary  conditions, 
where  either  u or  -gr  can  be  specified  on  any  edge  segment, 

for  regions  with  or’without  holes.  If  the  boundary  identity  (11) 
for  u is  expanded  in  a full  Fourier  series  in  the  range  (0,aq)  , eqn.(12) 
is  replaced  by  the  series 

* = I Bn  cosnx'  + I CN  sinnx^  = 0 ; n = 2irN/aq 

<1  n=o  H N=1 

Truncation  at  L harmonics  gives,  on  seroing  the  relevant  harmonics, 
BN  = o ; CN  = o ; N = 1 to  L 

Bq  = 0 ; (zero  harmonic  term) 

This  process  will  be  called  FuU_Harmonic  Matching  and 

, • 9T  + 1 pns  analogous  to  eqs.(16) , can  be 
the  corresponding  + i eqs.,  aiuit  & 


written  as 


Z At  r »•<*■>  -<«*;  * “k>S  - Sq(k;)cos(„x;«k)d^.o,  c 

i=l  0 


for  the  sequence  N“0  ; N = 1 with  k = 1 and  k , 

N =*  L with  k = 1 and  k = 2;  where 

a,  = 0 ; a0  = it/2 


» 


Each  side,  including  each  hole,  contributes  its  set  of  harmonic 
equations,  and  the  corresponding  cosine  and  sine  harmonics  from  all 
sides  are  grouped  together  in  harmonic  sets  N - 0,1 h,  to  constitute 

the  rows  of  the  coefficients  matrix.  3 

Again  as  in  (18)  convergence  of  series  09)  to  the  order  l/»  can 

be  insured  by  setting,  in  accordance  with  Fourier  series  theory, 

" ° l . 


= 0 and  x'  = a 


!!a  <*?.  o ( ' 

3x'  J 

A similar  analysis  applies  on  any  segment  q on  which^the  normal 
derivative  is  prescribed,  in  which  u is  replaced  by  « the 

boundary  identity  (10)  and  * (xj)  is  defined  accordingly. 

Four  vertex  equations  follow  from  eqs.(43)  for  each  boundary 
segment,  though  care  must  be  taken,  as  in  subroutine  COLMAT  in 
Appendix  B,  to  omit  equations  that  are  redundant  - as  is  the  case  when 

u is  continuous  at  a vertex.  _ # 

The  computation  of  eqs.(43)  for  all  possible  boundary  conditions 

requires  that  the  following  derivatives: 

3 3 . *2 

3yq  ’ 

be  available  for  all  functions  in  solution  (8) . 

We  now  proceed  to  reexamine  the  basic  problems,  and  to  obtain  the 

necessary  derivatives. 

gW 

(a]  Vertex  Functions  for  the  j vertex 

When  mixed  boundary  conditions  occur  on  the  sides  of  sector  iA 
Fig.3,  we  identify  four  possibilities,  distinguished  by  the  indicator 

NVEX  = NVER(j) , where 

(i)  Case  I:  NVER(j)  = 1 corresponding  to  conditions  u - 0 on  both 

sides  j and  (j-1)  of  sector,  which  as  in  (25)  gives 

A = o ; B - 1 ; \ = Wo  ; k = 

(it)  Case  XX:  NVER(j)  - 2 corresponding  to  conditions  n - 0 on  side 

j on  . 0,  and  |a  - o on  side  J-1  on  uhich  . 

the  normal  to  side  e - o and  can  be  denoted  by  the  unit  vector  . 


. * ’• 


Fig.  3 


Since  in  polar  coordinates 

it  follows  on  applying  to  eqn.(21)  that 

3u  = Xrx-1  [_a  sinA0  + B cosX0] 

3n 

Accordingly  on  analysing  the  sector  conditions,  as  in  (22). 
we  obtain 

A = 0 ; B = 1 ; cosXa  = 0. 

the  eigenvalues  of  which  are 

A = kit/ a - ir/2ct  ; k = 1,2,... 

ad  the  corresponding  corner  functions  are  given  by  for.  (21)  on 
using  the  above  values  for  A, B,  and  A.  _ 9u  __ 

(iii)  Case  III:  BVER(j)  - 3 corresponding  to  conditions  ,n  - 
6 . 0,  and  u - 0 on  0 - o.  Results  follow  similarly  and  are  give, 

in  table  I under. 


3u 

9n 


= 0 on 


(iv)  Case  IV:  NVER(j)  = 4 corresponding  to  conditions 

Q = 0,  and  ^ - 0 on  0 = a.  All  results  are  collected  together  xn 

table’ I und!?,  in  a form  that  facilitates  programming. 


NVEX 

Xk 

A 

B 

1 

ktf/a 

0 

1 

2 

kir/a  ~ w/  (2a) 

0 

1 

3 

kir/a  ~ it/ (2a) 

1 

0 

4 

kir/a 

1 

0 

(46) 

(47) 

(48) 

(49) 


TABLE  I - VERTEX  FUNCTIONS 


» < -i,- 


mm 


The  corresponding  Vertex  Functions  for  j vertex  are  denoted 


by  V,  . where 
3 kj 


V.  . H u = r [A  cosX  0 + B sinX  9]; 
kj  k k 


The  physics  of  the  problem  requires  that  X^  be  positive.  The 
normal  derivative  for  any  direction  y^  follows,  ns  in  (46): 

3u  A.  „ A,  r3u  A 1 3u  A, 

' Vu  ' ' % r * 7 ae  31 

q x. 

= {-Asin(A ' 0 + 6 . ) + B cos(X'0  +<j>  .)}Xr  » X1  = X 1 

qj  qj 


A -A . 

q j 


on  noting  that 


y'  . r = cos(ir/2  + ®)  = -sin(<fr^.-  9) 

O'  . $ = cos(<J>  . - 9) 

q qj 

where  y'  is  a unit  vector,  normal  to  direction  x',  which  makes  an 
angle  A with  y axis,  or  an  angle  ir/2  + <j>  . with  the  jth  side. 

q qj 

[b]  (Full)  Kdge-Functions 

Since  expansion  (39)  involves  a full  Fourier  Series  expansion 
on  the  side  AB  in  Fig.l  with  periodic  extension  to  all  other  points 
on  the  infinite  line  of  which  AB  is  part,  it  is  appropriate  to 
introduce  periodicity  into  boundary  conditions  (35)  for  the  half- 
strip problem  of  Fig. 2.  Accordingly  we  take  the  boundary  conditions 
on  the  infinite  sides  of  the  strip  as  , 

u(o,y!)  = u(aj *yj)  ; yj  > °> 

which  on  applying  to  solution  (35a)  require 
m = 2irM/aj 

Superposition  (37)  must  then  be  replaced  by 

OO  , 00  I 

u = A + I A .e_myj  cosmx!  + Z Be  ®yjsin  rax' 

0 M=1  mj  j M=1  mj  3 

where  Fourier  integrals,  similar  to  (38),  follow  for  A^  and  in 
an  actual  half-strip  problem.  As  previously  we  associate  the  cosine 
and  sine  edge-functions  in  (54)  with  the  corresponding  cosine  and 
sine  harmonics,  (M  corresponds  to  N)  in  expansion  (39)  - and  arrange 
that  the  corresponding  columns  and  rows  intersect  on  the  diagonal 
of  the  coefficients  matrix. 

The  above  cosine  and  sine  edge-functions  are  combined  in  the 
following  definition: 


, V \ •*  ■ • V , 


+J'i* 


( Fit  11)  Edge -Functi  ons  , for  side  j. 


E .=E1.+E2.  *,  Ek. 
mj  mj  mj  mj 


Ak.  c~mjyJ  cos(rn.x!+  a.);  k = 1,2 
mj  J J K 


2nM/ a j ; 


= 0 , <*2  - /2  > 


where  k = 1 and  k = 2 give  the  respective  cosine  and  sine  edge- 
functions.  The  terra  (Full)  is  usually  omitted,  as  expansion 
(39)  is  more  economical  computerwise  than  the  sine  series 
expansion  in  (12) . 

If  truncation  on  any  side  j in  expansion  (39)  is  taken 
at  L.  then  L (Full)  Edge-functions  from  side  j should^  included  in 
solution  "mix"  (8)  to  match  up  the  corresponding  harmonics  for 

N _ if2, L.  This  leaves  the  zero  harmonics  in  (41)  - the 

. total  number  for  the  problem  being  tallied  as  NZERO  - unmatched 
by  corresponding  edge-functions,  and  an  additional  NZERO  Vertex 
Functions  must  be  provided  to  match  these.  The  equations  for 
the  zero  harmonics  are  inserted,  as  rows  in  the  matrix,  directly 

after  harmonic  set  L. 

As  in  eqn. (51)  the  required  normal  derivative  for  the 
direction  y'  follows  as: 

q 


8Ek. 

"Li  - 


V(Ek.) 

mj 

3Ek. 

1 j3x'. 

J J 


'j 


= -Ak,  m.  e m3yj  cos(ra  x'.+  $ , + «k>  ; ♦qj“*q  *j  » 


since 


= "sin^qj  J yq  • yj  = C°S<i>qj' 


[c]  Curved  Boundaries  Including  Holes 

The  polar  solutions  of  Laplace's  eqn.  in  (21)  apply  to  closed 
curves  - where  0 acts  as  a parameter  for  points  in  the  boundary  with 
a range  (0,2n)  - provided  that  A is  an  integer,  say  A = k.  Fositive_ 
values  for  k' correspond  to  an  outer  boundary,  and  negative  values  to 

an  inner  one. 

Development  follows  as  above  for  (lull)  Edge— Functions , where  the 
range  of  2tr  for  0 corresponds  to  a^  for  x^ . Accordingly  when  truncation 
is  at  L,  for  any  closed  curve  we  must  include  in  set(8)  the  functions 


{ 


r { F.  p.  ;} 

j=i  k=o  J 

where  the  asterisk  denotes  that  the  terms  only  apply  to  closed  curved 
boundaries  j arising  in  the  problem,  with 

Pkj  = *\lj  cos(MjG)  + rMj  sin(Mj0) 

M = 1,2,3 .. 

Mj=  M,  if  j relates  to  a curved  outer  boundary 
=-M,  if  j relates  to  a curved  inner  boundary 
When  k = 0,  we  require  the  zero  order  harmonics 

Poj  = Aoj  + Boj  log  r 

where  the  log  r term  is  omitted  if  the  section  is  solid. 

The  origin  for  the  coordinates  (r,Q)  can  be  taken  at  any 
. suitable  point  within  the  corresponding  closed  curve,  and  could  be  taken 
at  different  locations  if  both  inner  and  outer  boundaries  are 
closed  curves. 


(57) 


(58) 


(59) 


The  normal  derivative 


follows  from  eqn.(51)  on  putting  X = N^. 


Note:  Functions  (58)  may  be  thought  of  as  Edge-Functions  for  the 

corresponding  curved  boundaries,  and  must  be  inserted  in  the 
appropriate  rows  in  the  harmonic  sets  in  the  Matrix. 


(d]  Logarithmic  Singularities  at  Vertices:  Log-Vertex  Functions  ( 

As  noted  in  the  analysis  of  vertex  functions  (21)  infinite 

derivatives  of  order  p,  where  p > X,  , occur  when  X,  is  not  an 

K K * * 
integer.  On  performing  a limiting  analysis  as  X^  ->  K , where  K 

is  an  integer,  it  can  be  shown  that  the  solution  to  eqn.(20) 

k k 

denoted  by  , for  X = X^  = K , where 

V*.  = -|r  [Ar*  cosXO  + Br^sinXe], 

KJ  o A 

= [Ar^  cosX0  + Br*  sinX0]log  r 

-AXr^sinXG  + BXr^cosX0, 
satisfies  the  following  conditions:- 

(i)  Function  (60)  satisfies  Laplace  eqn.(20)  since  the  derivative 
w.r.t.  the  parameter  X of  any  solution  (21)  of  eqn.(20)  also 
satisfies  the  same  equation. 


(60 

(61 


21. 


(ii)  the  coefficient  of  log  r - the  singular  part  of  the  solution  is 
ZCro  on  the  sides  0 - 0 and  0 - a,  since  A,  B and  X are  determined  in 
Table  I for  the  specified  boundary  conditions  on  the  angle.  However, 
in  the  interior  of  the  angle,  0 < 0 < a,  the  coefficient  of  log  r is 
not  zero , and  hence  the  solution  provides  logarithmic  singular  behaviour 

within  the  angle.  * 

(iii)  an  analysis,  as  in  (27),  of.  the  total  heat  flow  from  term  Vk] 

through  any  circular  arc  isolating  the  vertex,  gives  - on  interchanging 
the  operations  of  differentiation  and  integration  - 

q .ml-  (rx(l  - cosAci)  ] 

J 3a 

= rX[(l  " cos  Act)  log  r + AsinAa] , x = Xjc  '» 
and  consequently  form  (60)  is  physically  admissable  in  solution  W(8) 
provided  condition  (28),  A.  > 0,  is  satisfied. 

Note:  Functions  V*.  will  subsequently  be  called  LpG-VERTEX_FUNCTI^. 


[e]  Derivatives  required  in  set  (44) 


(62) 


The  dorivMiTT^r  follows  from  those  given  in  eels.  (51)  and  (56)  on 

replacing  * by  - W2,  since  the  direction  is  »/2  behind  that  of 

A, 

V 


On  applying  the  operator 


3 

3x' 


= xr  . V, 

q 


to  eqs . (51)  and  (56)  respectively  we  obtain 


2 • y 

9 Vkl  = {-A  sin(A"  0 + 24  .)+B  cos(A"0  + 2$  j)}A  X'r 
3x' 3y  V 

* q A"  = A~2 ; A ' =A-1 


2 


= Ak.  e“mjyi  sin(m.x!  + 24  . + «k) 

3x'3y'  mj  J J J 4J 


q q 


(63 


Derivatives  of  fnnetions  (58)  folio,  from  these  in  (63)  on  potting 

^Difficulties  are  encountered  in  deriving  analytical  expressions  for 
darivativos  Eor  Vertex  Bqs.C'O)  in  the  case  of  ourved  boundary  segments, 
as  shown  in  section  3(a),  and  numerical  differentiation  is  recoma, ended 

in  these  cases. 


i 


Accordingly  to  enable  the  program  LAPGEN  in  Appendix  2 to  deal 

with  all  boundary  segments,  whether  straight  or  curved,  numerical 

differentiation  is  used  for  derivatives  9'^q  required  in  Vertex  Eqs.(45) 

3x‘ 

q 

for  all  straight  line  segments.  This  is  based  on  the  two-point  formula 


f'(x) 


f (x  + Ax)  - f (x  - Ax) 


and  the  corresponding  computation  - involving  the  computation  of  two 
point  values  - we  shall  refer  to  as  a "two-point"  equation.  Similarly 
analytical  expressions  for  the  derivatives  of  V . in  (60)  are  best 
obtained  by  computing  numerically,  using  formula  (64),  the  of  the 
corresponding  derivatives  of  V,  ..  These  are  given  in  POLW,  where  the 
increment  AX  = delta  is  read  in. 


PROGRAM  l.APGEN 


The  program  in  Appendix  A(l)  now  requires  to  be  expanded  to  include: 

(i)  Examine  all  vertex  eqs.  for  redundancy,  and  - for  convenience  - 
set  up  all  information  to  enable  the  corresponding  point-eqs.  to 
be  put  up  in  rows  of  matrix,  after  those  for  the  zero  harmonics. 

The  assigned  subscripted  variables  are  - position,  (XFIX,YFIX); 
associated  slope,  BFIX;  side,  JFIX;  function  indicator  MBFIX; 
location  subscript,  (LOC) ; and  position  on  side,  (DF IX).  The 
number  of  equations  is  tallied  by  NFP. 

(ii)  Total  no.  Vertex  Functions  NORF  required  is 

NORF  = NFP  + NZERO  (66) 

where  NZERO  is  total  no.  zero  harmonics  (-  no.  of  sides  NS  in  a 
Laplace  problem) . 

(iii)  Provide  for  distribution  of  NORF  Vertex  Functions  between  the 

th 

vertices  of  polygon,  the  number  required  at  j vertex  being 
computed  as  NAD(j).  Locate  any  integer  values  for  A and 
arrange  to  include  the  corresponding  Log-Vertex  functions  (60) . 

(iv)  Set  up,  in  subscripted  locations,  the  i nf  o irmci  tier,  required  for 

computing  the  corresponding  column  sets,  NP,  of  the  matrix:- 

function  type  NFN;  vertex  location  HA;  associated  A and  B 

from  Table  1 in  coef(l,NP),  coef(2,NP);  and  eigenvalue  A ^ in 

E.  It  is  economical  to  put  up  functions  (52)  and  (57)  two 

columns  at  the  time  - corresponding  to  each  value  of  M - and 

each  two  such  columns  is  said  to  constitute  a colunm  set. 

' » 

Vertex  functions  (51)  are  put  up  in  single  columns,  or  colunm 
sets  of  one. 

(v)  Put  up  right  hand  side  of  matrix  as  its  last  column  set,  LP. 

The  above  tasks  are  performed  in  subroutine  COLMAT  in  Appendix  B(l). 

Function  indicators,  MT,  are  introduced,  in  accordance  with  the 
following  table: 


MT 

1 

2 

11 

12 

3 

Function 

u 

9u 

9u 

D2u 

3yq 

3x<j 

Sx'  9y ' 

q q 

Table  2:  Function  Indicators 


Note  The  harmonic  conjugate  v of  u is  given  by  MT  3 3. 


A more  elaborate  solution  program  SOLCOR  is  required,  which  uses 
Gaussian  elimination  on  the  harmonic  eqs.,  and  single  pivoting  to  locate 
the  maximum  element  in  the  corresponding  column  for  zero  harmonic  and 
point-eqs.  These  latter  equations,  since  they  are  matched  by  Vertex 
functions,  cannot  be  put  up  in  a diagonally  dominant  pattern  in  the 
matrix. 

SOLCOR  solves  coefficients  eqs. (17)  not  only  for  the  specified 
truncation  level  at  L harmonics,  but  can  also  obtain, as  a bye  product  - 
requiring  only  a few  percent  increase  in  the  overall  computer  time  - 
the  solutions  for  respective  truncations  at  L-l,  L-2  and  L-3  harmonics. 

The  indicator  LS  specifies  the  number  of  solutions  required 
where  LS  < 4 in  the  present  program.  Accordingly,  exact  solutions  are 
obtained  to  LS  mathematical  models,  of  increasing  complexity,  for  the 
given  problem. 

Provision  has  to  be  made  to  provide  the  derivatives  listed  in 
set  (44)  for  the  point-fns.  computed  in  subroutines  EDGEF  and  POLW. 

Subroutines  POLC  and  PMAP  are  added  to  provide  respective  Curved  Edge 
Functions  for  circular  or  elliptic  boundary  segments  and/or  corresponding 
shaped  holes. 

Consequential  adjustments  are  required  in  LAPEX  in  the 
information  input,  and  in  the  putting  up  of  the  matrix.  Points  on 
ellipses  are  located,  if  required,  by  the  usual  parametric  coordinates 
x = a cos  <j>,  y = b sin  <|>,  and  are  computed,  together  with  corresponding 
slopes,  in  subroutine  ELPS. 

t 

Provision  is  made  to  adjust  the  number  of  integration  points, 

NOK(j),  on  any  side  by  reading  in  a proportionality  factor  NHS(j). 

These  are  roughly  proportional  to  the  lengths  of  the  corresponding 
boundary  segments,  and 

NOK(j)  = 2*FSET*L*NHS ( j ) + 3 (67) 

where  FSF.T  is  a factor  > 1.  Fset  increases  proportionally  the  number 
of  integration  points  on  all  sides,  as  is  required  later  when  using 
Harmonic  Fitting,  as  developed  in  Section  5,  in  place  of  Harmonic 
Matching. 


An  indicator  NBDY  is  used  to  indicate  the  type  of  boundary 
conditions  (2)  arising  in  the  problem.  Thus  NBDY  = 1 indicates  that 
on  the  j*"''  side  - where  the  occurance  of  non-zero  boundary  conditions 
is  indicated  by  the  subscripted  variable  NBY(j)  = 1 - the  boundary 
conditions  are  of  the  polynomial  form 

6 k-1 

8 : (x>y)  = £ C(j,k)t  , (65) 

J k=l 

where  t = xj/aj.  The  required  coefficients  C(j  ,k)  are  read  in  for 
each  side  which  has  NBY(j)  = 1. 

Computation  of  point  values  of  (65)  is  done  in  a subroutine 
POLP,  which  also  provides  numerical  coefficients  - without  requiring 
to  read  them  in  as  in  the  case  of  NBDY  = 1 - for  the  torsion  problem  as 
indicated  by  NBDY  = 2.  Additional  cases  can  be  inserted,  corresponding 
to  values  of  NBDY  > 2,  in  POLP  as  the  user  requires.  • . 

The  checking  of  the  boundary  residuals  is  as  before  except  that 
it  is  done  for  LS  truncation  levels. 

Production  is  based  on  one  data  card  for  each  line  of  production, 
but  is  now  done  in  a special  subroutine  PRGDN , which  enables  various 
operations  indicated  by  NCODE  = 1,2,3  - differentiations,  principal 
stress  computations  and  integrations  - to  be  performed.  The  user  can 
readily  incorporate  into  PRODN  any  further  facilities  that  may  be 
required. 

The  expanded  program,  LAPGEN,  with  explanatory  captions  is 
attached  in  appendix  B(l),  and  an  illustrative  example  is  given  in 
Appendix  B(2).  A user  guide  to  inputting  data  on  the  black  box 
principle  is  given  in  Appendix  C. 
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Fig.  4 

(a)  Curved  Edge-Functions  for  Curved  Boundary  Segment  AB 

If  the  mapping 

C = F(z)  (6S) 

maps  the  rectilinear  region  A'B'r'S'  in  the  £-plane  onto  the  curvilinear 
region  ABRS  in  the  z~plane,  then  any  solution  of  Laplace's  eq.(7)  in  the 
t;-plane  satisfies  Laplace's  eq.  in  the  z-plane. 

Analogous  to  the  half-strip  approximation  (35)  to  basic  problem 

(b)  at  (20) , the  solution  near  A*B*  - where  the  other  sides  are  assumed 
to  have  little  effect  - has  the  characteristics  of  the  solution  in  ABRS 
in  the  region  A*B*,  when,  as  in  (52) , we  set 

u(S,  nL)  = u(5,  n2)  ; K > (69) 

as  the  conditions  on  sides  BR  and  AS  respectively. 

These  conditions  can  be  satisfied  by  requiring  that  u be  periodic  in  n, 
the  period  being  n2  “ 0.*  1°  the  transformed  problem  in  the 

>*J 

C-plane,  the  solution  in  the  region  of  A B is  characterised  by  edge-' 
functions,  analogous  to  functions  (55),  which  may  be  written  in  the  form 


r 


T 


12  k k ~mi^i 

C . = C .+C  C.  = A . e J J cos (ro.n . +ot.  ) ; k = 1,2 
m.l  nj  mj’  mj  inj  J J k ’ 

nij  = 2nM/(n2-n1)  \ = 0,  «2  = n/2 

til 

where  the  mapped  plane  C corresponding  to  the  j 1 side,  when  curved, 
is  denoted  by 

h - V ini- 

and  may  be  termed  the  Sj-plane. 

- In  subsequent  work  it  is  advantageous  to  work  with  the  complex 

variable  £ ^ , and  to  write  function  (70)  in  the  form 

-m.£ . 


-m.£ . 

C . = Re{A  . e J J} 
mj  mj 


(70) 


(71) 


J 


[A1.  cos(m.n  •)  “A".  sin(ra.Ti-)]  (72) 


J J raj 


J J 


where  A . is  a complex  constant 
mJ 

i ti 

A . = A . - iA  . 

mj  mj  mj 


(73) 


and  Re  denotes  the  real  part. 

Accordingly,  the  appropriate  functions  must  be  introduced  into 
solution  mix  (8)  for  each  curved  edge-segment. 

We  now  proceed  to  find  the  derivatives  of  as  required  in  set  (44). 

Analogous  to  mapping  (68)  we  can  write  the  mapping  for  curved 
side  j as 

cj  " Fj<z) 

Since  C . is  a function  of  z and  z,  say 
mj 


Cmj  = G(z>z)  ; z = x+iy  » 


z = x-iy 


we  require  the  operators 

3 ^ 

ax' 

q 


q iL  + e1<f>q  i_ 

32  3z 


(74) 


= ie 


% 3 . 3 

H -7T-  ~ ie  n — 

32  3z 


(75) 


The  operator  follows  on  noting,  from  fig.l,  that  the  coordinates 

q 

of  any  point  P(x,y)  are  given  by: 


X = X + 

X* 

cos  (!) 

- y’ 

sin 

q 

q 

q 

q 

q 

y = y + 

x' 

sin  (J) 

+ y’ 

cos  <J> 

(76) 

q 

q 

q 

q » 

• ...  ■* 


23. 


where  the  coordinates  of  the  q^  vertex  are  (x^,y^).  On  operating 
on  G(z,z)  it  follows  that 

3G  = 3G  3z  DC  3z 

3x'  3z  3x'  + - 3x'  ’ 

q q 3z  q 

and  on  substituting 

3*  = 3*  + i iz  = e^q  • ii.  = e~l4>q 

3x'  3x'  3x'  ’ 3x'  ’ 

q q q q 

3 3 

the  operator  for  — , is  obtained.  The  operator  -g— , follows  similarly 
Xq  yq 

Accordingly  we  obtain 


3C  . _ -m.£. 

. . El  = Re  [a  ~ (a  J J11 

3y ' KGlAmj  3y ' U jJ 

q q 

i<J>  dc.  -m . £ . 

= Re[-im.  A . e q — J-  e ^ 
J mj  dz 


The  parameter  q^,  for  segment  q,  now  corresponds  to  - the  parameter 


for  points  on  qth  line  segment  - as  used  in  harmonic  matching  eqs.(41). 
Analogous  to  the  vertex  eqs.(43),  we  require  derivatives  w.r.t.  q for 

q 

all  functions  in  solution  "mix"  and  for  all  their  normal  derivatives. 
These  are  difficult  to  compute  and  are  best  obtained  by  numerical 
differentiation  using  "two-point"  eqs.  formula  (64)  as  set  up  in  LAPGEN. 


(b)  Example;  Elliptic  indentation 


If  AB  in  Fig. 4 is  a segment  of  an  ellipse,  with  principal  axes 
x'  and  y’,  the  required  mapping  is 

z’  = c cosh  C;  z'  = x'  + iy' 
from  which  the  point  correspondences  are 
x'  = c cosh  £ cos  q 
y ' = c sinh  £ sin  q 

The  mapping  is  rendered  single  valued  by  taking  the  limits 
£(0,«>)  and  q(0,27r).The  corresponding  singular  points  in  the  z-plane 
are  at  F'(-c,0),  F(c,0)  with  a cut  from  -c  to  c on  the  x'  axis. 


(78) 


(79) 


(80) 


On  eliminating  n from  eqs(80)  it  follows  that 


2 ,2r 
c cosh  F, 


c sinh  F, 


1 , 


(81) 


or  any  line  l = l in  the  £-plane  corresponds  to  an  ellipse  in  the 
z* -plane  with  semi-axes  a = c cosh  b = sinh  c = a - b . 

Similarly  n - Uq  corresponds  to  a branch  of  a hyperbola  that  has 
the  same  focii  as  the  ellipse.  Then,  it  is  easily  shown  that  AERS 
is  mapped  onto  A'B'R'S'  in  a one-to-one  manner. 

i it 

To  obtain  point-values  for  the  coefficients  of  A . and  A . 

nij  mj 

from  cmj » or  its  derivatives,  we  require  to  find  £ and  n for  any 
point  P(x',y')  in  the  domain  -A-  - v/here  (x',y')  are  coordinates  of 
P w.r.t.  principal  axes  for  elliptic  indentation.  Equation  (81) 
gives  the  quadratic 

2„,2  . 2 ,2  , 2.  ,2  „ .,2 
cl  + T(c  -x'  -y’  )-y’  = 0 ; T = sinh  F,  , 


from  which 


, 2r  _ 2 .2  ,2.  J.  2 ,2  ,2,2  . . 2 ,2 

2c  f = -(c  -x  -y'  ) + y(c  -x  -y  ) + 4c  y , 


the  positive  square  root  being  taken  since  T > 0.  When  y'  = 0 then 

„ '2,2,  '2  2 

T = x /c  -1;  x > c 

- n ’2  ~ 2 

- 0 ; x < c 


It  follows  that 


£ = sinh  1 (/t)  = loge(/T  + /T+r) , 

and  eqs  (80)  then  give 

-1  ,y'cosh 

n = tan  \r~-, — ■■--.  -■■)  , • 

where  n must  be  adjusted,for  quadrant  location  of  P,  to  lie  in  (0,2tt) 

A suitable  subroutine,  PMAP,  similar  to  EDOEF  and  POLC  in  Appendix 
B,  is  included.  The  values  for  and  n2  are  obtained  from  (84). 


Note  that  singular  points  of  mappings  (74)  are  not  allowable  in  the 
domain  of  solution  (8),  which  must  be  regular  at  all  points  other 
than  the  singular  points  specified  as  in  (3) . If  the  singular  points  of 
any  mapping  occur  in  the  domain  SL  of  the  problem  then  either  (i)  split 
J'Linto  elements  to  ensure  that  the  domain  of  any  set  of  curved  edge- 
functions  does  not  include  the  corresponding  singular  points  of  the 
transformation, or  (ii)  superpose  other  suitable  functions,  as  over  in 
(e) , to  neutralise  any  unwanted  singularities. 


(82) 

(82a) 


(83) 

(84) 


(c)  Circular  Indentations 


This  follows  as  in  (b)  on  using  the  transformation 
z - e*’ 

Alternatively,  if  AB  in  Fig. 4 is  an  arc  of  a circle,  the  corresponding 
basic  problem  can  be  taken  as  that  for  a curvilinear  region  bounded  by 
radii  r = r^,  r = r and  radial  lines  0 = 0^,  0 = 0^,  where  -*■  » 

The  solution  then  follows  in  form  (57)  on  setting 


(85) 


M.  = -2uM/(0  -0  ) ' 


(86) 


(d)  Smooth  Inner  Boundaries 

If  the  inner  boundary  is  a smooth  curve  the  corresponding  Curved 
Edge  Functions,  analogous  to  polar  functions  (53),  follow  from  eq.(32), 
on  adding  the  periodicity  requirement 

= ° ; n2  = 2tt  (37) 


(e)  Solid  Elliptical  Section 

Let  us  examine  the  mapping 

z = c cosh  C,  (88) 

as  in  (b)^  in  the  case  of  a region  bounded  by  a solid  ellipse,  shown  in 
Fig.  5. 


The  corresponding  Curved  Edge-Functions  for  £ =*  £ analogous  to  those 
in  (72),  that  decay  inwards  from  the  boundary,  are 

CM=  ^ ^ = e^,C  cos  Mr1  - s in  MrJ 


(89) 


Y 


Those  require  a cut  FF'  in  the  z-plnne,  which  corresponds  to  the 
limiting  ellipse  £ = 0.  On  taking  E(0,“),  n(0,2ir)  to  render  mapping  single- 
valued n is  discontinuous  on  the  cut,  having  values,  at  respective 
neighbouring  points,  (i)  of  q = a and  n = 2t r— g:  at  and  P^  and  (ii)  of 
q = iT-a  and  n = if+ct  at  P^  and  P^  where  0 < a < ir/2. 

It  follows  that  cos  Mn  is  continuous  on  cut  FF'  but  sin  Mq  is 
discontinuous  as  it  changes  sign  - and  accordingly  C is  discontinuous. 

The  functions  got  on  replacing  M by  -M  in  (89) 


C-M  “ Re{A-M 


e = e M^[A'  cos  Mq  + A"  sin  Mq], 

"-M 


(90) 


have  a similar  discontinuity  on  FF',  and  on  superposing, the  resulting 
discontinuity  on  £ = 0 is 

2 (A”  - Aiy  sin  Ha, 

and  this  reduces  to  zero  on  setting 

Analogous  to  eq.(78),  and  noting  from  eq.(83)  that 

Jc  = 

dz  c sinh  £ ’ 

the  discontinuity  between  points  P and  P in  the  derivation  of 

Cfrl  + C-M  is'  on  5 = 0:  q 


A" 

-M 


i(J> 

Re[iMe  q ' 


. 

ic  srn 


q>  {AM 


~ A_m  e"^}] 


(91) 


(92) 


2M  ,A  _A  ^ ,cos_Mn  + +A  ^ sinjln^ 


Re[e  4 <Va-m>  tStt  + (Va-m> 


sin  n 


(93) 


and  hence  the  discontinuity  in  the  term  cos  Mq/sin  q vanishes  if 


-M 


Also  note  that  at  the  singular  points  of  mapping  (88)  the  function 

g 

Cw  + C and  its  derivative  — — f-  is  everyiAere  finite  since  limit 
M “M  dy 

'q 

of  (sin  Mq/sin  q)  ■+■  M at  points  F’,  q = ir,  and  F,  q = 0 or  2u. 

Similarly  all  the  higher  derivatives  of  C +C  ,,  can  be  shown 
• M -M 

to  be  finite  and  continuous  at  all  points  within  and  on  the  elliptical 
boundary. 

If  a region  is  bounded  internally  and  externally  by  ellipses, 
then  the  curved  Edge  Functions  for  the  outer  bounded  must  be  rendered 
finite  and  continuous  by  superposing  the  corresponding  C when  the 
focii  F and  F'  of  the  outer  ellipse  are  in  the  region  of  the  problem. 


• N I ^ A .»* 


4 


(94) 


(f)  Approximate  Curved  Edge-Functions 

. . * 

If  the  mapping  (79)  is  used  to  map  the  side  A’B  of  the 

rectangle  A'B'R'S'  onto  an  elliptic  arc  in  the  z-plane  that  is  not 

too  different  from  the  curved  boundary  segment  AB,  then  the  functions 

C . as  got  from  the  elliptic  mapping  approximates  pretty  well  the 
^ _ -k  iV 

characteristics  of  the  solution  in  region  A B and  can  be  included 

in  solution  "mix"  (8). 

Foints  on  the  actual  boundary  can  be  expressed  in  terms  of 
£.  and  p,  where  n is  in  The  parameter  n can  be  used  for  the 

actual  boundary  points,  the  corresponding  £ being  where  the  cur/e  p 
intersects  the  actual  boundary  - and  p , for  segment  q,  then  corresponds 

to  parameter  x'  in  eqs.(51)  for  harmonic  matching. 

9 

Accordingly  on  exact  mapping  (79)  is  not  essential,  and  in  most 
practical  cases  a "fitted"  elliptic  mapping  is  sufficient. 


In  practice  this  would  require  the  fitting,  by  a least  squares  criterion, 
of  an  elliptic  arc  - involving  five  parameters  (axes,  centre,  orientation)  - 


to  the  prescribed  arc  AB.  We  shall  refer  to  the  resulting  mapping  as  a 
"fitted"  mapping. 


SECTION  4 - OTHER  BASIC  PROBLEMS 


(a)  RE-ENTRANT  ANGLE 


The  semi-finite,  or  half-strip  model,  used  in  arriving  at  Edge- 
Functions  (37)  or  (55)  does  not  apply  to  the  sides  of  a reentrant  angle, 
as  at  D in'Fig.6.  On  considering  the  side  DE,  the  assumption  that  the 
solution  on  DE'  is  related  to  a periodic  repetition  of  that  on  DE  is  * 
clearly  untenable.  Also  in  the  region  yj  < 0 the  exponential  becomes 
positive  and  obviously  cannot  represent  a propogation  into  the  interior 
of  the  boundary  actions  on  DE. 

However  if  we  divide  up  Fig.  6 into  two  convex  regions,  I and  II, 
by  taking  a cut  from  D anywhere  within  the  angel  E'D  C' , say  along  DA, 
this  will  ensure  that  y^  > 0 for  each  region,  and  that  all  the  corresponding 
Edge-Functions  have  the  characteristic  negative  exponentials  at  all  points 
in  their  respective  regions. 


34. 


The  solution  "mix"  for  eq.(l)  can  Chen  bo  written  as 

u = u + u , . . 

o (1) 

= u + u , . . . 
o (n) 

for  regions  I and  II  respectively,  where  denotes  all  the  functions 

that  are  common  to  both  regions  and  are  continuous,  and  with 

continuous  derivatives,  across  the  cut. 

The  functions  u,..,  and  u....,  which  are  confined  to  regions  I 

(i)  (n) 

and  II  respectively,  must  then  be  matched  across  the  cut  to  give 

continuity,  in  a least  squares  sense,  across  the  cut  for  u and  — . 

on 

Since  the  Edge-Functions  from  DC  and  DE  produce  a discontinuity 

across  AD,  this  must  be  negatived  by  superposing  functions  that  are 

characteristic  of  a line.  All,  of  discontinuities  in  u and  — . Analogous 
3n  o 

to  the  half-strip  model  (52)  the  solution  characteristics  in  region 
A*D*  are  similar  to  those  for  an  infinite  strip  bounded  by  infinite 
sides  perpendicular  to  AD  through  A and  D respectively. 

On  taking  edge-axes  x'  and  y'  for  AD  as  shown  in  Fig. 6,  we 
require  that  the  solution  u(x',y')  satisfy  the  following  boundary 
conditions  for  the  strip: 

(i)  u(x',y')  + 0 as  y'  •>  “ and  y’  -+  -®, 

(ii)  u(0,y')  = u(a',yr),  for  all  points  on  x’=  0 and  x’  = a1,  where 
a'  = AD. 

Since  u(x',y')  satisfies  Laplace's  eqn.  suitable  forms  for  u,  • 
in  semi-finite  strips  y'  > 0 and  y'  < 0,  analogous  to  forms  (54),  are, 
with  m = 2irM/a'  : 


u = E e 


-my 


{A  cos  mx'  + B sin  mx'} 


y'  > 0 


E e 


my 


{A1  cos  mx’  + B'  sin  mx'} 
m m 


< 0 


In  an  actual  strip  problem,  if  the  discontinuities  in  u and 
specified  as  f(x')  and  g(x')  respectively,  it  follows  that 


du 

3y' 


were 


(93) 


(94) 


f(x')  = E{  (A  -A')cos  mx'  + (B  -B')sin  mx'} 

in  m mm 

g(x')  = E{-m(A  iA')c os  mx'  + -m(B  +B')  sin  mx'}, 

° mm  mm  * 

and  on  using  Fourier  coefficients  formulae  the  coefficients  in  (94) 

follow. 

However  in  an  actual  problem,  as  in  Fig.  6,  the  coefficients  must 
be  determined  as  in  ( 11)  which  allows  for  the  influence  of  all  sides,  and 
functions  (94)  must  be  included  in  solution  "mix"  (8). 
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These  latter  can  be  identified  as  Edge-Functions  of 
type  (55)  for  the  upper  and  lower  side  of  the  cut  AD,  where  the  upper  side 
has  indicator  j and  the  lower  has  ( j + 1) . 


(b)  SINGULARITY  ON  EDGE 

Consider  a singularity,  as  in  (3), 

u = Q <5(n-n  ) 
o 

acting  at  the  point  P,  parameter  q , on  curved  boundary  segment  AB  in 
Fig.  4.  On  taking  ABRS  as  the  basic  region  with  conditions  (69),  and 
superposing  solutions  (72),  to  satisfy  the  conditions  on  AB,  £ = 
we  require 

/1\  (J\ 

Q 6(n~n  ) = I e [A' . cos(mn)  + A . sin(mn)  ] •• 

o mj  mj 

The  coefficients  follow,  on  taking  a summation  M = 0 to  M = 00 , as 


A<!>  e"”'1 

mj 

(1) 


n2~r)l 


/ 5(n-u  )cos(mn)dn. 


O 

(2) 


which  gives  for  Av . and,  similarly,  for  A 

mj  me  ml 

A(l)  _ Q e icos(mno) 
mj 


a(2)  = Q e 
mj 


n2_nl 
m£i 


sin(mr|  ) 
o 


n2-nl 


(95) 


(96) 


Series  (95)  provides  a useful  particular  integral  (4)  for  above  type 
of  boundary  singularity.  Since  its  value  is  known, Q5 (n-nQ) , on  AB, 
the  series  only  requires  to  be  evaluated  at  interior  points  (C,n) 
where  £ > and  it  is  then  dominated  by  the  negative  exponential 

factor 

Note  (i)  that  all  quantities  ,n^  >n2’n0  an<*  m re^ate  to  si-de  j» 
but  subscript  j is  omitted  here,  as  not  necessary. 

(ii)  If  the  singularity  acts  on  a straight  side  we  require  the 
obvious  replacements 

n *»  x' 


q = 0;  = 0; 


n2  = al5 


no  = a’ 


the  singularity  being  at  the  point  = a. 


*-  » 
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SUCTION  5 - DISCRETE  OR'l'HOCONAI,  FUNCTIONS 


If  the  function  y=f(x)  is  given  as  the  set  of  n discrete  points  x^ 

in  the  interval  (a,b) , by 

yx  = f(x^)  ; A = 1 ,n 

and  {0^(x))is  a set  of  m discrete  orthogonal  functions,  then  we  can 

approximate  f(x)  by  the  polynomial 

ra 

y = 1 Ck0k(x) 

A=1  k K 

by  minimising  the  sum  of  the  squares,  ql 2,  of  the  weighted  errors 
weight  factors  w^  - at  the  n data  points  x^: 

2 n - 2 
0 = X=!  WA(yA"yA) 

n m 2 

= E WAtyA_  1 Ck0k(XA)] 

A-l  A A k=l  k k A 

where  y,  is  the  value  of  y at  x . 

J A A 


On  minimising  q2  with  respect  to  the  unknowns  c^;j  l,...,m  by 


setting 


we  obtain 


fl-Oi 

oc . 

J 


0 ; j = 1,. . . ,m;  n > m ; 


Z wA[yX  " E CA(X  )]  W = ° 

A=1  X X k-1  k 3 

If  the  set  {0  (x)}  is  orthogonal  over  the  discrete  points  with 
k 

weight  factors  w^,  this  means  that 

Z wi»k(xl)0J(x1)  - »j  4 ; 

A J. 

where  N.  is  the  'norm'  of  the  function  0 j , and  6^  i.s  the  Kronecker 
delta.  Accordingly,  equation  (A-5)  gives  on  replacing  subscript  j 


by  k: 


l E WAyA0k(xA) 

k a=1 


r 
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If  the  set  0 (x)  is  complex,  then  0.(x)  in  (101)  and  (102)  must  be 
• I 

replaced  by  its  complex  conjugate  0 ^ (x)  since  Nj  must  be  real  giving 


E,  wA(xx)0j(V  = Vk 

X=1  J 


1 n * 

Ck  = N.  * VxW’ 

k A=1 


(103) 


on  noting  that  the  right  hand  side  of  (98)  must  be  real  which  implies 

2 

that  c^0^(x)  must  be  real,  and  that  n is  real. 

Two  limiting  cases  for  n are  of  special  interest: 

2 

(i)  If  n=m,  then  the  obvious  minimum  condition  for  n is  that  each  of 

the  errors  (y  -y  ) should  be  zero,  or  the  approximating  curve  y should 
A A 

pass  through  the  set  of  points  (x  ,y  ) , X=l,...,n,  and  hence  y and  y 

A A 

would  be  point-matched,  or  collocated,  at  the  above  set  of  points.  If 
a unique  solution  exists  it  is  necessary  that  n > m. 

(ii)  If  n ->  « and  m ->  then  the  discrete  set  x,  becomes  the  continuous 

A 

x in  the  interval  (a,b) , and  (98)  becomes  an  expansion  in  the  infinite 
set  of  orthogonal  functions  (0^(x)}.  When  0^(x)  is  complex,  the 
coefficients  follow  from  (103)  as 


Ck  ~ N, 


/ w(x)y(x)0*(x)dx, 


(104) 


and  the  condition  for  orthogonality  of  the  function  0^(x)  in  (a,b) 
follows  from  (102)  as  ^ 

Nk6k  = J w(x)0k(x)0- 

a J 

where  w(x)  is  the  associated  weight  factor  in  (a,b) . 


(105) 


Discrete  Fourier  Polynomials 
ikx 

If  0k(x)  = e in  (98)  it  is  easily  shown  that  these  functions 
satisfy  the  orthogonality  relations 


\ eikx  e"ijkdx  - 2x6 j ; i = /-l, 


(106)  . jj 


the  quadrature  expression  for  which  is 


• a 


T,t  I w.elkxX  e_ljXX  = 2rr6:)  , 


ivw  X=-n 


(107) 


* 
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where  the  weight  factors  are  w^ , the  origin  for  x being  taken  at  the 
mid-point  of  the  interval  (— ir , tt ) , which  is  divided  for  convenience 
into  2n  equal  parts. 

lkxX 

Result  (107)  shows  that  the  functions  e are  orthogonal  over 
the  infinite  set  of  discrete  points  x in  (-it, a)  with  weight  factors 

A 

w^ . We  now  ask  can  some  w^  and  distribution  of  points  x^  be  found, 
resulting  from  the  application  of  some  quadrature  formula  to  (106) , 
which  will  satisfy  relation  (107)  for  finite  values  of  n and  thereby 
provide  a method  of  deducing  a set  of  discrete  orthogonal  polynomials 
from  a corresponding  set  of  continuous  orthogonal  functions. 

On  investigating  the  trapezoidal  formula,  with  its  associated 
weight  factors 

w.  = j ; X = + n 
= 1 ; I A I < n 

and  taking  the  equidistant  distribution  for  x 


x.  = — 
A n 


; A = -n  ( 1 ) n , 


we  are  led  to  examine  the  finite  analogue,  S,  of  series  (107): 

n • /.  , v n . , A 

Q = V „ l(k“j)x.  = V „ „lX0 


where 


9 = it (k-j) / n . 


On  substituting  into  (109)  for  w^  from  (108)  and  applying  the  formula 
for  summing  (2n+l)  terms  of  a geometrical  progression,  it  follows, 

when  k i j,  that  i(2n+l)9 

S = -l(e~ln°  + ein°)  + e_ln°  , 


where 


inO  ivr(k-j)  . . . , . ... 

e = e = cos  (k-j)  it  = + 1 ; k £ j , 


depending  on  whether  (k-j)  is  even  or  odd.  Accordingly  it  follows 
that  S = 0. 

If  k = j,  then  the  series  (109)  reduces  to 
n 

S = I w = 2n, 

•v  A 


and  hence  for  all  k and  j 


=>  2n  6:' 


(112) 


V 
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which  establishes  condition  (103)  with  N,  = 2n  for  the  functions  *• 

ikx  k 

e A when  k is  any  integer,  including  zero.  It  follows  from  (103) 

that 

, n -ikx 

^ r-  A 

ci  _ E e 

k 2n  , \ 

A=-n 

and  if  we  take  a correspondingly  balanced  polynomial  form  for  (98) 
by  writing 

m 

^ ikx 

y - E c « , 

k=-m 

it  follows  from  (113)  that 

* 

C-k  = Ck  » 

A 

where  c^  in  the  complex  conjugate  of  c^,  and  on  writing 

c.  = a.  + ib. 
k k k 

series  (114)  reduces  to  the  real  series 

m m 

y = c + 2X  a,  cos  kx  - 2X  b,  sin  kx, 

° k-i  k k-i  k 

where  on  taking  real  and  imaginary  parts  of  result  (113),  we  obtain 

1 n 

2ak  = n , Z WAyACOS  kxA 
A=-n 

1 n 

2b,  X w, y,  sin  kx, 

k n , A A A 

A=-n 

On  setting  x = x'  - tr,  we  easily  show  that  (115)  and  (116) 
are  invariant  in  x’.  Hence  the  above  formulae,  hold  for  the  x’  interval 
(0,',tr),  and  this  can  be  mapped  onto  the  interval  (0,c)  by  setting  * 

Xs  = 2ttx'/c,  replacing  kx'  by  2irkx"/c  Or,  on  dropping  the  dashes,  this 
means  that  (114)  and  (116)  apply  to  the  interval  (0,c)  on  replacing  k 
by  2irk/c. 


Trigonometric  Interpolation  Series 

Since  boundary  identity  (39)  can  be  approximate  by  the 
truncated  series 

Li  Li 

tfi(x')  X cos  nx'  + X C„  sin  nx'  , 
q N =0  N q N-l  N q 

in  the  interval,  0 < x'  < a ; n = 2rrN/a  , we  might  regard  the 

- q - q q 

trigonometric  series  on  the  r.h.s.  as  a trigonometric  interpolation 

scries  for  ^'(x')  in  the  specified  interval.  The  fitting  coefficients 
q 

tlien  fol low  by  discrete  least  squares  as  in  eq.(98),  the  coefficients 
being  given  by  formulae  (116). 


We  observe  that  the  series  (115)  is  similar  to  a Fourier  series, 
truncated  at  m terms,  for  the  interval  (-ir,n)  and  that  the  summation 
formulae  (116)  for  the  Fourier  coefficients  2a^  and  -2bj_  are  what 
would  be  obtained  by  evaluating  the  corresponding  Fourier  series 
integrals,  analogous  to  these  in  (16),  by  the  trapezoidal  rule, 
involving  the  division  of  the  interval  (~7r,7r)  into  2n  equal  parts. 
Hence  a working  rule  for  the  harmonic  equations  in  sets  (16)  or  (41) 
Evaluate  all  Fourier  integrals  by  the  trapezoidal  rule 

This,  in  effect,  substitutes  series  (116)  for  the  corresponding 
integrals,  as  is  required  when  harmonic  matching  is  rep1 aced  by 
discrete  least  squares  minimisation  of  the  boundary  residuals  on 
each  boundary  segment. 


vi 


f 
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SECTION  a - DISTINCTIVE  FEATURES 


The  distinctive  features  of  the  Edge -I-’unction  Method,  as  illustrated 
by  the  examples  in  Appendices  A and  B and  the  formulation  in  the  present 
paper,  are: 

(a)  Algebra : 

The  functions  introduced  in  solution  "mix"  (8)  may  appear  complicated 
at  first  sight,  but  an  effective  algorithmic  method  of  controlling  the 
resulting  algebra  - called  the  Computer  Form  Method  - was  developed  by 
Quinlan  [12,13].  Accordingly,  Edge-Functions,  Vertex  Functions  and  any  of 
their  derived  functions  like  normal  slopes,  moments  or  shears,  can,  on 
calling  the  appropriate  subroutine,  be  obtained  as  readily  as  any  trigono- 
metric or  exponential  function. 

(b)  Programming : 

A systems  approach  to  programming,  based  on  a main  program  Quinp 
together  with  18  subroutines  - each  with  its  own  definite  task  to  perform  - 
is  given  in  [4].  This  has  been  used  to  advantage  by  subsequent  research 
workers,  Tai,  Mash,  Dashmukh,  O'Callaghan  and  others,  to  considerably 
simplify  and  shorten  the  programming  tasks  arising  in  extending  the  Edge- 
Function  Method  to  problems  of  greater  complexity,  in  vibrations  and 
shallow  shells.  The  program  Quinp  [4]  and  its  accompanying  program- 
description  is  still  adequate,  though  it  is  proposed  to  issue  shortly  an 
up-dated  version  to  include  cuts,  cracks  and  reentrant  angles. 

The  program  LAPGEN  given  in  Appendix  B(l)  illustrates  the  chief 

# . . » . * 
features  of  QUINP.  If  one  is  familiar  with  LAPGEN,  the  various  parts  of 

QUINP  should  then  be  readily  understandable. 

(c)  Computing-Time : 

Considerable  computing  time  is  saved  by  using  the  discrete  - 
rather  than  the  continuous  - least  squares  method  of  minimising  the 
boundary  residuals,  as  developed  in  Section  5,  thus  saving  considerable 
time  for  setting  up  the  coefficients  matrix.  The  rows  arc  arranged  in 
the  coefficients  matrix  so  that  the  equations  for  each  harmonic  are 
grouped  together  in  successive  bands,  thereby  producing  a strongly 
d Lagonal i sed  system  with  a cons iu  bie  saving  in  solution  time.  More- 
over this  arrangement  of  the  matrix  enables  the  solutions  corresponding 
to  a lesser  number  of  harmonics  to  be  deduced  without  any  appreciable 
increase  in  computer  time.  Accordingly,  as  a routine,  solution  vectors 
corresponding  to  terminating  the  boundary  identities  at  L-3,  L-2,  L-l  and 
L hnrmoni.es  respectively  are  computed  and  tested  for  each  problem. 


-/  — -»■  — . _ «■ . — r- 
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(d)  Acco.ptabi  1 igy  of  Solution: 

Each  solution  vector  provides  an  "exact"  solution  to  a problem 
governed  by  the  same  differential  equation  but  with  slightly  different 
boundary  conditions  to  those  specified.  Such  solutions  may  be  regarded 
as  mathematical  models  of  the  physical  problem.  In  each  case  the 
difference  - termed  the  Boundary  Residuals  - between  the  boundary  values 
as  computed  and  the  specified  values,  is  computed  and  reported  through 
its  approximate  root  mean  square  value  - r.m.s.  - on  each  side  of  the 
boundary  and  for  each  boundary  condition.  The  set  of  r.m.s.  values 
provides  a simple  yet  comprehensive  reliability  test  to  enable  an 
engineer  to  decide  whether,  or  not,  to  accept  the  results. 


(e)  Convergence  Demonstration: 

A practical  "convergence"  demonstration  is  provided  by  the  routine 
provision  in  each  problem  of  several  solutions,  corresponding  to  increasing- 
harmonics  and  matrix  size,  with  r.m.s.  values  presented  for  the  resulting 
boundary  residuals.  These  invariable  decrease  rapidly  as  the  number  of 


harmonics  used  increased.  Likewise  product 


result^  f 


(c.g.  normal  slopes, 


shears,  harmonic  conjugate,  etc.)  as  required  at  interior  points,  are 
alx/ays  computed  for  a number  of  different  harmonics  and  hence  as  in  examples 
given  in  appendix  B their  "convergence"  can  be  seen  at  a glance.  Mo  other 
competing  system  - finite  element,  finite  differences  or  boundary  integral  - 
can  offer  this  effective  comparison  of  the  effect  on  the  results  of 
increased  computer  time  expenditure,  without  involving  a very  considerable 
increase  in  computer  time  over  that  which  would  be  required  for  a one-shot 
solution  based  on  the  largest  matrix  size  involved. 


(f )  No  pre-computer  processing  of  the  problem  is  required 

As  can  be  seen  from  Appendix  C and  from  the  data  cards  at  the  end 
of  each  example,  no  pre-computer  processing  is  required.  Only  the  geometrical 
and  load  data  and  material's  moduli  are  required,  together  with  one  control 
card.  These  do  not  require  any  knowledge  of  E.F.M.  for  their  preparation. 
Accordingly  E.F.M.  can  be  operated  completely  as  a "Black  Box".  This  is  in 
sharp  contrast  to  the  intricate  element  networks  required  in  finite  element 
and  boundary  integral  methods.  Production  is  based  on  a single  data  card 
for  each  production  set,  requiring  the  computation  of  a specified  function 
at  a specified  number  of  equidistant  points  on  a specified  line. 
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APPENDIX  A - LAPEX 


Program  LAPEX  consists  of  a simple  main  program  LAPEX  with  four 
subroutines  EDGE!’,  QPOLAR,  POI.W  and  QSOLVE(NE).  The  program  is 
described  briefly  at  the  end  of  Section  1 and  is  attached  together  with 
illustrative  example  (37)  as  Appendices  A(l)  and  A(2) . 

The  corresponding  data  cards  are 

(1)  CONTROL  CARD 

L,  Ns,  Ndiv,  Mdiv,  Nprog 
where 

'J 

•'  L = Truncation  level 

Ns  = No  sides 

Ndiv  = No. divisions  used  in  integral  evaluation 
Mdiv  = No. of  check  points  on  sides 
Nprog  = Example  no. 

(2)  COORDINATES  VERTICES  OF  POLYGON 
x(j),  y (j) 

(3)  COEFFICIENTS  FOR  BOUNDARY  CONDITIONS  (37)  on  sides 

c(j,k)  (k  = 1,4)  - for  each  side 

« 

(4)  PRODUCTION 

ui.  V1#  U2,  V2,  Mdiv, 

for  function  u at  Mdiv  equidistant  points  on  line  (u  ,V^)  to  (U2,V2) 

It  is  urged  that  the  reader  should  become  fully  conversant  with 
LAPEX,  before  proceeding  to  its  fuller  development  in  LAPGEN  in  Appendix 


45. 


(5)  Two  Minor  Subroutines  are  included:  QPOLAR,  to  determine  polar 

coordinates,  and  CAUSS  to  determine  division  points  and  weight 
factors  for  both  Harmonic  Fitting  and  Gaussian  Integration.  The 
program  LAPGEN  follows  with  numerous  explanatory  captions,  the 
references  being  to  the  main  paper.  . 

An  illustrative  example  is  appended,  as  Appendix  B(2),  to  attached 
program  LAPGEN,  Appendix  B(l) . This  deals  with  the  torsion  of  a quadilateral 
se.ction  with  an  elliptical  cavity,  Fig.  7,  for  four  different  truncation 
levels . 

An  effort  has  been  made  to  make  the  output  in  Appendix  B(2)  self 
explanatory.  The  Boundary  residuals  are  less  than  1%,  and  the  differences 
between  computed  quantities  for  levels  LL  2 and  LL  = 4 are  seen  to  be 
considerably  less.  ■'  * 

Solutions  of  acceptable  engineering  accuracy  are  provided  by 
truncation  level  LL  = 1 , involving  only  53  equations  with  a time  requirement 
on  any  computer  of  less  than  150%  of  the  time  it  would  require  to  solve  53  !' 

linear  equations  using  Gaussian  elimination. 
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APPENDIX  !J  - I.APGEN 


This  program  consists  of  main  program  LAPCEN  and  10  subroutines 


consisting  of 


COLMAT  Sets  up  data  for  points  at  vertex  j in  the  vertex  equations 
taking  care  to  eliminate  any  redundant  equations.  It  then  assigns 
the  functions  for  solution  "mix"  (8) , and  sets  up  the  necessary 
data  for  their  evaluation  at  any  point.  Assigned  points  and 
functions  are  then  printed  out  as  a useful  aid  in  checking. 

POINT  FUNCTIONS  for  indicators  MT  = 1,2,3 

1 2 

(a)  EDGEF  Evaluates  coeffs  of  unknowns  E . and  E . for  Edge- 

mj  mj 

Functions  in  eqs  (55)  and  (56).  . 

(b)  POLW  Evaluates  Vercex  Function  V,  . in  eq  (50)  and  also  Log- 

. * 

Vertex  Functions  ^ , eq  (60),  where  they  replace  . 

(c)  POLC  Evaluates  Harmonic  Polars,  eq  (58),  including  the  zero 
polar  log  r. 

(d)  PMAP  Evaluates  curved  Edge-Functions  based  on  section  3 for 
mapping  z = c cosh  w as  given  by  eq  (89)  including  the  special  case 
of  a solid  ellipse. 

(e)  POPP  Evaluates  r.h.s.  of  matrix.  Particular  integrals  not 
yet  included.  Provision  is  made  for  boundary  conditions  of  types 
NBDY  = 1 and  NBDY  =2. 

S0LC0R  Solver  routine  for  LS  truncation  levels,  based  on  single 
pivoting  for  columns  using  the  relevant  equations. 


PRO DM  Arranges  for  production,  other  than  computation  of  results 
for  MT  = 1,2,3  as  indicated  by  Ncode  = 1.  This  is  divided  into 
two  categories:  Ncode  =2,  and  Ncode  = 3 corresponding  to  differen- 

tiation and  integration  respectively  as  the  main  operations  involved. 
The  several  cases  in  each  category  are  designated  by  the  indicator 
NOFN.  Provision  is  made  for  torsional  rigidity,  shear  stress, 

resultant  stress  and  shear  lines  for  the  torsional  problem  developed 
in  Sokolnikoff  [14].  Other  cases  can  be  added  as  they  arise,  and 
thus  a comprehensive  program  for  Lnplacc/Poisson  problems  can  bo  built 
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APPENDIX  C - DATA  INPUT  FOR  LAPCEN 


The  data  sets  required  are 


(1)  CONTROL  CARD 

L,  NS,  NB,  NPROG,  NBDY,  LS,  NMAT,  NPRIN,  MDIV,  DELTA,  FSET , CMOD 
where 


L 


NS 

NB 

NPROG  = 
NBDY  => 


LS 

NMAT  -> 
NPRIN  =* 
MDIV  => 
DELTA  = 
FSET 
GMOD 


Maximum  no.  harmonics  - if  set  as  percentage  of  optimum, 
program  will  set  corresponding  L. 

No.  boundary  segments 
No.  closed  boundaries 
No.  assigned  to  problem 

Type  indicator  for  non-zero  boundary  conditions 

6 

1 indicates  polynomial  Z C(j,k)t 

k=l 

2 indicates  torsion  conditions  as  set  in  POL?. 

No.  Comparative  solutions  required 

Matrix  print  out  given  if  set  NMAT  = 1 

Cives  print  out  of  residuals  on  boundary  if  set  NPRIN  = 1 
No.  of  checking  points  required  on  boundary 
Increment  for  use  in  numerical  differentiation 
Degree  of  smoothing  required 

Torsional  rigidity;  if  dealing  with  torsion  problem 
read  in  value  for  GMOD. 


Note  L,  LS,  MDIV,  DELTA  and  FSET:  if  any  of  these  are  left  blank,  program 

will  provide  an  appropriate  value.  Consequently  only  NS  and  NB  must  be 
specif ied . 

(2)  DATA  FOR  EACH  SEGMENT 

MB(J),  NHS(J),  NBY(J) , NTYP(J) , X(J) , Y(J) 
where 

MB(J)  = Function  indicator  MT  as  defined  in  Table  2 
NHS(J)  = Proportionality  factor  to  regulate  no.  fitting  points 
per  side,  approximately  proportional  to  length.  If 
left  blank,  will  be  set  to  one. 

MBY(J)  ^ Indicates  if  set  to  one  that  non-zoro  boundary  conditions 
<y[  type  NBDY  - 1 occur  on  side  j,  and  corresponding  coeff 
must  be  read  in  as  data  set  (4)  under. 
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NTYP(J)  => 


X(J),Y(J)  => 


Boundary  curve  type  indicator  for  segment  J where: 

0 indicates  straight  line 

1 indicates  curved  indentation 

2 indicates  curved  mound 

3 indicates  closed  curve 

Coordinates  of  vertex  J,  where  appropriate 


(3)  ADDITIONAL  DATA  CARD  FOR  EACH  CURVED  SEGMENT 

AE(  J)  , BE  (J)  , XE  ( J)  , YE  ( J)  , 7.E(J), 
where  for  elliptic,  or  circular  boundaries 

Semi-axes  are  AE,  BE;  centre  (XE ,YE)  and  inclination  of  major 
axis  ZE  (as  multiple  of  ir/2), 

(4)  NON-ZERO  BOUNDARY  COEFFS.  FOR  NBDY  = 1 

C(J,K),  (K  - 1,6)  , for  sides  J with  indicator  NBY(J)  = 1. 


(5)  PRODUCTION  - one  card  per  production  set,  unless  when  NCODE  = 2 

or  NCODE  = 3 a second  card  must  follow  to  describe  the  function 
that  is  sought. 

U , Vx,  U2,  V , BXK,  MDIV,  MT,  J,  NCODE,  NORN,  NAXIS, 
where  production  is  for 

(i)  SEGMENT  (U.^V  ) to  (U2,V2)  or  if  J is  specified  it  is  for 

segment  J of  boundary 

(ii)  BXK  = angle  associated  with  production  (i.e.  shear  • 

axis)  as  a multiple  of  7r/2, 

(iii)  MDIV  = No.  equidistant  points  at  which  results  are 

required . 

(iv)  NCODE  =>  Indicates  type  of  production  required 

= 1 : Functions  MT  = 1,2,3  as  in  Table  2 

= 2 : Functions  involving  some  differentiations 

- 3 : Functions  involving  integration 

(v)  NOFN  •-->  Indicates  different  functions  available  for  each 

of  the  production  types  NCODE  = 2 and  NCODE  = 3. 
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6-  07 

14 

- .9  24l-  ■ 

0 . 1 1 4 7 L — p 3 

i c 

- .2  2VA--  • 

• 1.72.  - - > 

16 

. i : 

- . . 7 .? 3E-  . 2 

1 7 

- .13  7 : —OF 

. 1 ,74.  - .6 

17 

- . 1 3 24  -3 

- 0 . 6 7 5 

1 7 

- . 1 .<•  7 . - - • . 

3.1  1 4 L - 

2C 

.2133 

- . 136-02 

? 1 

- .1311- 

. 1 J9L--6 

2? 

- . 5 6 p - 

- /.  5m  1 66  - -2 

2 3 

- .247  . - ■ 

. . ' - -7, 

if  : if  if  - if  if  if  if  if  * * ‘ 
■:•  : ir- if  =1=  i-  * ♦ * * * 

J 

3 

It  - > 5 .1 

. ' I : ) 

Z '_■<()  It-  ",  D I V 

3 

1 062  .-06 

3 172-2'? 

2 . 33C-26 

4 c 92 u -03 
3: 63l-: 6 
1 3 *:  - 3 

3 1 601:  -06  2 

, 1 4 3 . - C’  3 - 0 

, 2 3596-98 
17 1 -04 
,2  32  7 -06 
,47176-04 
7 - 

. 3 1 7.  - 4 
, 3 7 74  _ -to 
.27420-04  -0 

. MnZ  6 

.274  -04  - 


$ « $ * * «**  : - <•  ■=  *##* 

^ , C 1 I JM  u ’ j I M T 3 u J 


. 74432 
. 1 )2 7. 
.3771. 
.14  4 

>•  $ <:  <•  $ if 

IV:  ' J L 


-0  6 0 
-05 
-0  6 

-04  — 0 

-06 

ijt  ^ s'*,  jjc  * 

I ML 


6 1 E - 0 5 

L FI  t L A 30 
4 

121  -n6 
, r 76 6L-'''  3 
, 1 7666-^6 
, 1 4 3 2 _ - ' 3 
l o r-rc 

7:,/'  -C 

,24146-04 
,2 1622 -0  3 
. ; 1 • -,'7 

.3  1':.-"; 

-04 

. 

, 3 35  7 0 -"4 

. 

.6463  - 
. 1 74  -07 

. 3 ; > ~ " 4 

’ - ' 

.127 

.447  -• 7 

.5134  - 
. 4 O . -97 


,•>  * j -,jT  ' 0 4 - A U i . V 1 t U2  * VZ  . 

. ' 3 . ' J j ' 0 ■'  J 0 . 1 " 0 ..  1 

if  $ fc  i-  ].  • I-  >-  -1-  *>  A — ^ fc  >:  ■ ’>  : — # — — if  if  £ if  i 

'TS  ON  Ll  J 1 * V l 

1 . 4 0 v - 7 


f '4  * 1 


•i 


1* 


* * It  * £*>!:*  t * * * *#  * ■ <:  .:  t(t  if  t<  ft  ft.  ft  ft.  ft  <1  * ft  ft  ft  ft  ft  ft.  ft  ft  ft  * ft  ft  ft  ft  ft  ft  ft  ft  ft  * 


1 


3 

A 

5 

6 

7 

3 

3 
13 
1 1 
12 
1 1 
1 A 


C LAP3EN  PRDG  .AM  FjR  ifci-4-  .:  AL  LAPLACE  PROBLEMS  APPLMDIX  Ml* 

ftftftft.ftftftftft.-ftftftftft.ft  ftft.ftftftft'-ftftftftftftftftftiOeftftftif'ftiftft'cftftftiCcftftftftftitcft#** 

C 3 *.  A 3 M X ( 1 i ) , V ( 1 1)  , C ( ) , 6 ) , A ( U)  > :'  ( 2 » 1 C ) , A N G ( 1 0 ) iLI  1 3 D ) , ■■'  K S ( A , 1C 

CG(  9 9,  10  D > , ML  l A ,UK  ) , A-  (ll),Btlll),Xr  ( l 1 ) , Y ( 1 1 ) . XF  I X { 5E  ) * YF  I X ( 5 

C,XCL( 2,  -0  ) , YEL  ( 2,  53*  , . M?  , 9-  ) , >AR  ( 2 , 1 I , 0 v CD , DF  I X ( 5 0 ) , 1 1 1R  ( A * 

C'JMMJN  'F  I <(  3C  ) ,ZL;  ( 1 1 ) , M A ) ,C  Jr  F ( 2, 1 3 } , M , DF  L T A , X C AR  ( AC  ) , A 5j  ) 

ZD „ ( 1 DfNAI  1 i , . , 

C " ".Ia  ; ( HI  fMH»NA  (11).  .J.  ( 1 i ) ,NEQ  < A)  , LP , NAX  IS»  NEL  l* 

: MMON  N ] K ( 1 1 ) , I*  ( 1 1 ; . . H , ( 1 L ) , . Y!  11  ) , T Y P 1 1 1 ) t JF I X { 5 ),M  F I X ( 5(J 

ZD'*  MM  NFP  ,NZlR  1,  MS,  , NULL ( 1 . » ,L  C<5.  I « *JBDY  t NF  POL 
[ A*  ' 1 3 1 DM  OLD  ( AJ  * A 1 ) 

31 ''.MSI  Dm  j A DO  Y ( l 1A92)  , I JK(  A 70) 
l J J W A _ M C . ( P A . . D Y ( 1 ) » X ( 1 ) ) »(  I J K C 1)  i iF  L (1)  ) 

J 1 2 .3  A I - , 1 1 A 32 
° A^DY ( 1 ) = " . 3 L + - . 

1 2 3 A CDNTINJL 

i O 12  36  1 - 1 » A I . 

12  3 5 I J < ( I ) = 0 


PRIM 


DAI  .. 


15 
1 6 

1 7 
18 
1 ) 

2 0 
21 
22 
2 3 
2 A 
23 
26 

2 7 
2 2. 

29 

30 

31 

3? 

33 

3 A 

3 5 
36 
3 7 
38 
3 7 
AO 
A 1 


A3 
AA 
A 5 
A 6 


MIT.  (6,  799) 

A R I T E ( S , 7 9 i ) 

CUNT  RCL  CAM 

• ••  (5,330  )L,NS,NB,NPROG,N  Y»  LS»  1MAT,NPR  IN.MDI  V,DELTA,rS  r ,3MC 

MITE  (6,79  90 ) NPROG 
H I T r ( t , 799) 

< LA  ( 5 , 7 -1  9 1 ) ( XC  AR  ( I ) ,1  = 12*9 
, RI ! L ( 6, A9  31  ) ( XCARl  I ) , I =1 , A ) 
ft  < I T - ( 6 , 7 9 3 ) 

, ’ IT; ( 6,5  1.  1 ) 

WRITE  ( 5 » 50  0 ) L , NS  » MB  , NP  OG , NB  D Y » L S , V AT  »NPRI M,  MO IV,  Di  LTA»FSE  I ,G* 

? 1=3. 1A159265 

ft  R I T C ( ‘ , 7 9 3 ) 

ft  T I HI  6,  A 9 91  ) NS 

X RITE(6*50_2) 

'101X05  = 3 
jn  l j = l , \ s 

5 u f NAT  A 

R L 1 p ( 5 t 799  9)  ( J ) , JH  3 ( J ) » N D Y ( J ) , NTYP(J)  , X ( J ) , Y ( J ) 

■'j:ar3s  = m:a  ’ j s + i 

NCAR1T  T LLIlS  THL  DAT  a CARDS  AS  E AD  IN. 

1 Ift  ; I T-.  ( 6,  A 9 )1 ) J,  •:  ( J ) , ,US(  J ) , >JPY  ( J ) ,NTYiM  J)  ,X  ( J ) , Y ( J ) 

"M  d 78  J = 1 . N S 

I F NHS(  J > L FT  F LANK  S-  T (v  1 

I F ( \ 3 S ( J * ) 7 6,877,378 
° 7 7 jHS ( J ) = 1 
376  CDNTINJE 
37  9 ,.R  I T-  ( 6, 79  3) 
i .D  IC  = 0 

3.  1233  J=;,NS 
I F ( NTYJ( J ) ) 1233, 1233,1232 

AJLIIUNAL  DATA  F 0 ' CJKYLD  SIDES 

123.  RCA!  (6,79?  ) A L ( J ) , b L ( J ) , XL  ( J ) , YE  ( J ) , ZL  ( J ) 

N C A 9 D S = 4C  A DS+  i 
1 F ( INDIO  ) U32,  1 ft  3 2 , 1231 
12  32  ftRI  TM  6,  7 9 <) 
ft  MTC(  5,A9  J2> 


4 8 

r.  3 IT:  (6,61  3) 

49 

12  3 1 

wRITE(6,4998)6L(J)tBL  ( J ) , X 2 ( J ) ,YE(J)  . Z T ( J ) 

SO 

Z E ( J ) = 3 . 6 *Z  r ( J ) 

.1 

1233 

:ont i n j e 

c 

1F  jrUA  4 ( < .1  left  lam-  sft  tj  j _>  j /» L VALjES 

62 

I F ( JELT Ml  37,  1 ,36,  12  i7 

5 1 

12  36 

J L 1 A = 0.  10002-03 

64 

1237 

I F ( 8 St  r 112  )9,  1238, 123  9 

55 

1238 

F SC  1 = 1.0  J 2 

66 

12  39 

IK  J 3 ) 1 240,  1240 ,124  1 

s r 

i 2 4 0 

lu  = 1 

66 

124  1 

i'jdi:  = o 

3 9 

I FM3DY-2  > 1245,  1244,  124  4 

60 

1244 

JO  1246  J = 1 , N S 

61 

124  6 

MdY(  J ) = i 

62 

03  13  1247 

63 

1245 

3 3 2 J = 1 , N S 

64 

I F M 3 Y { J ) ) ,2,3 

C 

.<  /\ , :u:  = FS.  C(  J,K)  , (K=l  ,5)  FUR  NB  PY  = 1 

66 

3 

RE  A (6,  30  Cl)  ( C ( J , K ) , K =■  1 * 6 ) 

66 

j:ar3s  = 'i:a  js+i 

6 7 

IF  ( MCI  C)  1 243,1  242,  1243 

68 

1242 

i md i : = i 

69 

R I r E { 3 , 7 9 i ) 

70 

/,  R I T b ( 6 , 4 9 ■ j 3 ) 

71 

124  3 

f.  \ I TE(  6, 5 0 1 ) J , ( S ( J , K ) , K = 1 ,6) 

72 

COMTIMJE 

73 

f,  I T E ( 6 , 7 9 )) 

r 

= in.  SHARP  CMNERS  30  3 0JN3A  RY 

r 

' ,UK  ( J ) = :JG.  I .TFGRATICN  n T S . 

c 

JRDWCJI  - Vi  HE  I HARMS.  FOR  SEGMENT  J BEGIN  I EACH 

c 

,V2RlJ)  = VERTEX  vD.  FOR  v.RTEX  J 

74 

12  4 7 

MSS  = N!S 

76 

I F(  V TY>  ( V S )-  3) 5 ,4,4 

7 fo 

4 

NSSMS-1 

7 7 

1 F f OS-2 ) 6, ; 0, 5 

78 

10 

N S S=  3 

7 9 

5 

•1R3W  ( 1 ) = 1 

80 

NH=  0 

3 1 

30  9 J = 1 , ':S 

32 

N 0 < ( J ) = 2*LMHS(  J>  *F  St  1 + 3-9 I Y P ( J ) / 3 

83 

M V7F  R ( J ) =3 

84 

MR-V(  J + l ) = RJVK  J ) + 2 * • IS  ( J ) 

85 

^HsNlH*2*£>n  ^.(J) 

36 

i r < \ r y 3 ( j ) -3 ) 6 , ,,9 

87 

6 

j j = j - 1 

38 

IF( J-l) 7, 7,8 

8 9 

7 

J J = MSS 

90 

8 

NVER  ( J ) = 2*(MB(  J )-  1 ) +MI  ( JJ  ) 

91 

9 

CON  T IN Jr 

c 

o 

SECTION  2.  SIDE  LENS  TH  S , A ( J)  , SL  OPES , B ( 1 , J 1 

c 

VERTEX  ANGLES, ANG ( J) ,F0R  POLYGON. 

ENSURE  ANMJ)  10  INTERVAL  C,2*Pl> 

IT  MSS)  16 , 16, 1 1 

1 1 xnss  + i ) = xt  i ) 

Y ( MSS  + 1 ) = Y(  1 ) 

J"1  1 2 J = 1 , MS 

; 1 ( 1 : r y >’  ( J ) ) 1 !> l . , l 61 u , 1 6; c 


4t&*' 


r ' — — * 


9 7 

1520 

CALL  t-PS ( J»0, J> 

93 

j3  TO  12 

9 9 

1510 

x a = x r j ♦ i ) - , ( j > 

13  J 

Y A = Y ( J ♦ l ) - Y ( J ) 

131 

: O3  3L  A ( X A * Y A , < R , 1-  A , P I ) 

132 

31  1 , J ) = B A 

133 

3 ( 2 * J ) = 3 A 

104 

A ( J ) =3  3 

13  5 

12 

continue 

106 

0 3 1363  J = 1 » MS 5 

10/ 

IF1 J-l ) 15, 13, 14 

IDS 

l 3 

AND  11)  = 3 I - ( l , 1 > + B(  2 , .SO  ) 

139 

jU  1 3 1 5 

113 

14 

A 5 3 ( J ) - 3 I - . ( 1 , J ) +B(  2 , J-  1 ) 

1 1 1 

1 5 

1 F ( A N J { J ) - *3  I ) 1 5 5 0 , 1 55  : , l 540 

112 

154 

A M.~  ( j ) = ANG  ( J ) - 2 *P  I 

11  3 

1550 

IF( ANSI  DI1  555,  1560 , 2 o 

114 

155  5 

ANSI  J > = A MCI  J ) +2 -PI 

115 

156  : 

E <3  N T 1 '4  J 

c 

AO  JUST  L 3 L S TO  APPRDRIAI 

L VALUE 

IF  NECES 

S ARY 

116 

1 F 1 L - 2 6)18,1/,!/ 

1 1 7 

17 

LMA  x = 1 l OO-  ,rP)  / ; S 

113 

.=( _1AX*,_ )/100 

1 1 9 

18 

I F ( L - 2 ) 19 ,20,20 

123 

1 9 

1 — *' 
L.  “ C 

121 

20 

IF! . -L5 )2  1 ,2500,2500 

122 

2 1 

lS=. 

123 

2 5 00 

IF(LS)?  2,22,23 

124 

22 

L S=  c 

125 

JO  13  2 3 

126 

2 3 

1 r ( L 5 - 4 ) 2 5 1 , 2 5 C 1 , 2 4 

12  / 

2 4 

L S=  4 

128 

250  1 

:dnti mj: 

C 

CALL  SUE  GUTI4L  C L MA T TO  SET  UP 

r i xeo  p 

or  NTS 

C 

amu  as s i y-i  r no  cdrrespc  . 

niNG  COLUMNS  OF 

’ AT  "I  I X 

129 

1 6 

CALL  COLMAi 

c 

SCC  T I 00  SLT  UP 

C CLUM';S 

ND  OF  ■••.A 

T K 1 x ; 

L 

CM  CYCL 

1=2  AND 

HARMON  I 

CS  ON 

r 

EVALJATE  INTEGRALS  CY  TRAPEZOIDAL  ■ 

. U L 5 WITH 

NUK ( J ) 

PT  S . P 

130 

M .,305  = : 

131 

JCC  = 0 

132 

N3P  = { N-I*L  ) /? 

1 33 

DO  51  \l3  = I .LP 

134 

4.  S = 1 

135 

11  ( N 3 - N P 3 ) ; 5 , 2 5 , 2 6 

136 

25 

MCOLS  > 

c 

: u L U 'I ;N  ..  AT  A AS  S.T  U9  I 

N COL  MAT 

1 3 / 

26 

J E = N A ( N P ) 

138 

NJF*9  = M=  J(  4T») 

139 

O’ni'Ui 

1 4 C 

X JE  =X ( J t ) 

141 

Y Jr  = Y 1 J 5 ) 

142 

2 J C = 3 ( l,JCi 

143 

JRO  = NH*  L 

144 

JRPr JR3+9Z  30 

1 

00  501  1=1,2 

146 

J J=  \iFP 

14/ 

mm:^=o 

143 

I F ( 1-1)28,  7,  28 

1 4 ) 
is: 

151 

is  2 

is  i 
1 s<, 
iss 

156 

157 
1 5 8 
157 
160 
161 
162 
163 
1 OS 

165 

166 

16  7 


27  JJ=\S 

2 H DO  5 00  J = 1 . J J 
If  ( I - 1)  U , ;D, 31 
3''  < < = VI3<  ( J ) 

GAL.  35.  JSStKK,NU'3L  ) 

YT=Y3(  J ) 

J P-  J 

OH  13  327 

31  K K.  = 1 

y=ldg(  j ) 

M T = V R F I < ( v ) 

I F(  Mr-  10)  327,32  ,,323 
72  8 K < = ? 

'U=^T-1D 

3 2 7 D3  53  < = 1 ,KK 

IF ( 1-2) 32, ^3,33 

I j TE  GR A T I 0 4 PL) I ; T S AND  ,._lGlT  FACTORS 

32  r<=ri<> 

„ K = a ( K ) 

J \l  3 = M T Y P ( J ) 


1 6 8 
167 

170 

171 

172 

173 
17  4 
175 
1 7fc 
177 
1 78 


1 79 

180 

1 B 1 

1R2 

18  3 

184 
135 
1 86 
137 

16  8 

IP  5 
170 

17  1 
1 9? 
1 7 3 
1 74 
173 

176 

177 

178 

1 7 y 

2 00 
201 


DATA  FOR  POINTS  ON  .LLIPTIC  Cb^VE  \'C.  NNCL 

IF { JM3  ) 330  , 23^ , 631 
33  1 IF ( < - 1 ) 332, 3 32  , 333 
33?  I F ( N5-l  ) 3 7 i,  334 ,33  1 
334  NNrL=NMEL+l 
7 Jr _ ( J ) =M Ni  L 
GALl  E-PSU, NINEL, KK) 

3 3 3 NEL=N3EL<  J ) 

XK  = X E L I NE  L , K ) 

Y<  = YEl.(  NEL,K  I 
3X<=3EL(NEL,K) 
j3  TO  34 
r 

DATA  F.:R  POImTS  i.M  .TRAIGHT  LIN.  SEGMENTS. 

33  0 X<  = X(J)+OK.v(X(J+l  ) — X ( J ) ) 

Y<  = Y(J)OK*(Y(J4l)-Y(J)) 

B X < = 3 { 1 , J ) 

GO  TO  34 

DATA  Fl  7 FIXED  POI.TS  AG  ASS1 GML . IN  COLMAT. 

3 3 -A  = _ 3 0 ( J )+K-  I 
X<  = Xr I X ( M ) 

Y < = Y f I < ( '3  ) 

J J = J F I X ( M ) 

B X<  - 3 F I X ( ) 

D<  = P- I X ( Y ) 


GALL  3UGP0UTI  IE  S FOR  POINT  VALUES 


34 

GO 

10 

( 35 

, >5  , 

36  , 

37,37 

, 1 7 0 ) , NFN 

3 5 

CA. 

. _ 

EOOE 

r ( xk 

, YK 

, B X K , 

XJ_  ,YJE  ,6 JE,CP,PI ,MT,VA,YB) 

GO 

TO 

38 

76 

GA.L 

POL5 

( <K  , 

YK  , 

B XX  , i 

K , FT,  .COL  E , J!’  , V A ) 

NCO-S 

= i 

J. 

GO 

TO 

33 

3 7 

C A = 

GO 

cF  ( 1 

, .P) 

Ct>  = G0EF(2,  P) 

C . L L J0L«  ( XK  , YK  ,B  XK  , XJL  , YJ  , 13  JE  , EM,  C A , CB  , NT  , 
G3  TO  33 
3 7 0 X J E 2 X E ( JE  ) 

Y J r = YE  ( JE  ) 
bj:  = 7f(  j:  i 


V A * P I , DEL  T A) 


202 

1FINFN- 6 >3/1,372,  3 12 

22  J 

3 72 

CAL.  P'3A3(X<,Y<,t3XK.,XJ.  ,YjL,t3J  . E M,  MT  , JC  , \l  A,  V H ) 

204 

J.l  T 3 38 

2 05 

37  1 

CAL.  PJuC  ( .•  < , YK  , XK  , •;  J . ,4  J ,B  Jb  , EM*  MT,  VA  * VB  » Pit 

22  b 

38 

FI  1 )=VA 

221 

r ( ? ) = vj 

c 

COMPUIC  WAF  lx  ENTRIES  A i J r I , S £ ’J  7 IJ  / ( J7 

, J < ) 

20b 

I F ( 1 - 1 ) 4 1 , 4 1 , 3 7 

c 

P 0 I , v T i < •[  ^ P(  i j|  • t .0  . 

209 

3 7 

J R = J <P+  J 

PIC 

J 2 43  J3~  i»vJ3Jl . 

PI  1 

jx=  j ; : i j: 

P 1 2 

4 0 

U(JR,JX)  = GU:<,JX)-F(JC)*1(-1)**K) 

211 

jO  T 3 53 

214 

4 1 

J < = J 3 J ♦ J 

c 

/ERL  4 A KMC  MC  E .3. 

215 

Da  42  j:=i,\ools 

P 1 6 

jx=jcc*  iz 

? 1 7 

4 2 

j(J,*JX)  = G(J  1,  J x)  + ,mk*F  ( JC) 

Pie: 

M=  Cr 

?l ) 

YM= MH  S ( J ) 

c 

C 3S  I ANU  SICE  ( IARMCJN  I C LOS. 

22  0 

DJ  44  M S= l * L 

221 

2 3 4 4 '3  = 1,.'  K 

222 

V1='NN  + 1 

223 

33  4 4 '12  3 = 1,2 

224 

[ .TA=(  \l  LG-  1 )*PI  *0. 51  x 

225 

J R = M S-  1 ) *.mH+  NR  GW  I J ) + ( 2*  I K-  i ) ) + NEC- 1 

226 

C:  = C3S(  2*  P I *JK  * iN  + eCTA  ) 

22  7 

D3  44  j;=i,n;ols 

22  o 

j x= j:c+  j: 

222 

44 

u(JR*JX)=G(JR»J<)  + rtK  *CG  -'F  l JC  ) 

222 

5 2 

; 33 l I NJ5 

231 

500 

comtimj: 

232 

50  1 

com  n nje 

233 

5 1 

jcc=jc:+mc  lc 

? 3 4 

NE'  = \’E*-1 

c 

P IK  I . I GUT  OF  CU2FFS.  MATRIX  IF  '.MAT  SET 

TO  1 

235 

IF ( M/A  1 I 3 42, 342 ,341 

23t» 

341 

« < 11  1 I 6 , 7 9/) 

237 

« R I T I (6,50 1 8 ) 

238 

33  347  j;=l,MEL 

23  7 

aR  I T:  ( 6,5016) JC 

240 

WR  I TC  ( 6, 50  ; 7 ) 1 ^ ( JR  , JC)  , JR  = 1 , M ) 

241 

3 4 7 

CQNT I N JC 

242 

a'  <1  T f ( 6 , 7 7 ) ) 

243 

342 

IF  ( NC-4J)  343 , 34  j,,  351 

244 

34  3 

D 3 3 4 b’  J 3 = 1 » N C 

244 

33  348  )C=.,ML. 

246 

348 

HDL  3 ( J K , JC ) =G ( J R, JC ) 

24  7 

a R I T (6,79/) 

c 

CALL  SUBROUTINE  TO  SOL  VE  MATRIX  FJR  LS 

soljt 

c 

a:ju  CHECK  '.CSULTS  if  : . E < 4 1 

245 

351 

cal.  S3_:n; 

2 4 7 

J ri  MC-40)  3 2,  352,  Jr>3 

? 30 

352 

a R I T C ( 5 , 7 9 ) ) 

251 

WRIT;(6,5"19)LS 

25  2 

a R I T _ I 6 » 5 C 4)(LL,LL  = 1,LS) 

25  j 

LJ  3500  J R = 1 * N E 

254 

33  350  LL=i,LS 

r U.)=C.  ' + 09 

:m  34  i j:  = 1, 

34  > KLL'^-ILLIHUL:  I J iJJ  *S3L  ( LL , JC  ) 
3 5 1 C DM  T I 'J  J 
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3D8 

3 F I X ( < ) =3  XF 

3D  i 

I'D  ? 1 9 J 5 = 1 , L P 

UD 

J ; = N M N* ) 

3 I 1 

NFN=NKJ(NP) 
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3?  1 

-J  ID  217 

322 

214 

^ A l _ P J i_  "^  ( ^ t ;<K  y UK  f If  f ? J ) 

32  3 

N ^ D l S = 1 

32^ 

33  ID  217 

32  5 

2 15 

C<=-DtF{],  . P ) 

326 
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cjnolate  coeff-s.  in  col.  jx 

c 

SDLlTIOM  SL)  L ( LL»  LX)  A . 4 D F l \ 

f 

SID  J AND  TRUNCATION  LEVEL 

3 3 8 

DC  218  J C = , 1COLS 

339 

JX= JCC+  JC 

34C 

D3  213  LL=1,LS 

341 

2 1 
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l F ( N C 3 3 C ) 2 2 1 9 , 2 2 1 7 , 2 2 1 

3*»<t 

2219 

DO  220  L L = 1 , L S 

34  4 

22  r 
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IE  .CODE  >1  FuA  1 1 r R PKOCtf  SSI 

c 

34  7 
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376 

N C A R 3 S = nG  A , 2 S+  1 

37  7 

/,  R I T r ( t , 4 9 y 1 ) ( XL  AO  ( I ) ,1=1,2. 
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L X<  = 5 X A *7  I : . 5r + 2 0 

38  7 

, s 1 3 : = j 

388 

IS  J=  1 

"589 

NC  Y=  1 
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c 

1JF  . = . . 
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52  8 

IF  ( < - 2 ) 3 5 , .5,54 

52  9 

P 4 

0<= 1 - D< 

530 

8 5 

orrx(M)  = n; 

531 

JF I X ( NX  ) = JJ 

53? 

Ms  c I X < N X ) = 3 ( J J ) + 1 

533 

I F ( J T Y3 ( J J ) ) 86, 26,3  7 

534 

8 6 

xnx(\3)  = x(jj)  + iK*(x( 

535 

YF  I X < NX  ) = Y<  J J ) + r.K  *(  Y( 

536 

FlX(NX)=E’ll,JJ) 

Q 

SIMILAR 

537 

TO  1 ID 

5 33 

8 7 

L 1 = <>AX  ( 1 , J J ) +DK* ( PaR  ( 

5 39 

LL-  L (01) 

540 

X l = :OG(  ' 1 ) 

541 

Yl  * SI N ( 71  ) 

542 

1HY1)  . ),  8 , 2 9 

543 

88 

Y 1 = .10-10 

544 

8 9 

25  = --t;  ( IJ  1 * X l/(  AE  < J J ) 

545 

4 n = 3l 

545 

1 F ( r ! - » I ) 9 1 , 9 1 , 9C 

54  7 

92 

A 9 0 = ? * P I 

54  3 

91 

1 * . 

54  9 

\ ^ r r> 

55C 

■Ml);-.'),  , < ,9 

i pt  i c s eg*  r s . 


iURRGUT 1 


r—- 


65  I 

92 

AD 9=  * I 

552 

9 i 

3 : ! X l N3  ) = 2 + 11  + AO 0 

5 5 I 

Xl  = Xl*4L'(JJ) 

56*, 

Y 1 = Y 1 «=  ( JJ  ) 

555 

XI  I X ( NR  ) = X (JJMX1*U.,(ZZ)-Y1*S 

556 

Yl  I XI  3 7 )-  Y ( J J)  +XL*S1  AM)  *Y1*  . 

55  7 

100 

contimj: 

558 

*> 

C-  u 

..UNI  I MJE 

55  9 

23  1 

contimj: 

C 

<Z  “i  V.  I ,Db  Ti.  ' I ,AC  Y I N 

c 

FONT T IOM  TO  ZCMC  AT  LAST 

c 

Mr  L I M _ o.  A , LlT T BY  COL. 

569 

10  1 

N L L I '1=1 

561 

00  1C3  3=1,  MS 

562 

I F ( Y i ( j ) - ? ( 1 0 2 , 10  3,  L 3 3 

J>5  J 

102 

\ L IM=J 

5 64 

10  3 

C 0 N T I M J r 

565 

IF  ( Mr  LI  A)  I 6,  10  5,  104 

566 

i :4 

M 1 = M < + 1 

56  7 

M l = M : +1 

568 

l o:  < m:  ) = mr 

56  9 

X F I X ( M 1 ) = X ( 2 ) 

6 7u 

V F I X ( N 9 ) = Y ( 2 ) ♦ B t ( 1 ) 

671 

Y F l X ( M X)  = I 

5 7 2 

Of  I X ( N3  > = 0 

573 

J F t X ( N < ) = 1 

c 

rt  mjvC  INDilT  R Y I M A C Y 1 

, T i J M U . » SI  =0  AT  T:  ;< 

574 

105 

, •=  M 7+  1 

575 

m:  = m : + 1 

576 

L J 0 ( M 2 ) = N R 

57  7 

jr [<i M3 1=  i 

576 

OF  IX ( M3)  = 0 

5 7 3 

MbF I X ( MR ) = 3 

59  0 

XF I XI  M3  ) = X ( 1 ) +A  n 1 ) 

53  1 

Y F I X ( M 3 ) = Y ( 1 ) 

592 

210 

Mr°=M: 

583 

J f 

5 84 

MZE30=MS 

535 

ME=ME+NZERi 

536 

M09F  = ME-J'H  3F 

p 

9U\I  = TOTAL  ML) . OF 

1Z_  0 = ( T M MG.  1 

53  7 

1 1 (NT  Y > ( 1 ) ) 23, 2 i,  2 1 

5Bt 

21 

NP=  y|»  + l 

I 3 COL  j.  FOR  CCNSTA.NT 

c 

ELLIPSE  3 LUO  FONCTIUN  FOl 

c 

I 1U1CA  in  BY  S: 

539 

NA  ( M •'’  ) = I 

5 3 

MFJ ( M5  3 = 5 

69  1 

..  ( M ’ ) = 0 

6 3 2 

MOn  = MO<F-l 

5 93 

in  MS-2  >24.23,2  3 

594 

2 3 

DO  240  3=  2 , IS 

596 

I T ( MTYn  J ) >240,240,23 

596 

? 3 ' 

J P=  M "*  ♦ 1 

5 3 7 

N A ( N 9 > = J 

5 3 tj 

. F L3  ( M 0 ) =5 

69  3 

if  ( A*s  ( J >-E,  ( J >-..  IOjL -04  >2  32 ,23; 

JtUM  4 . , PRO1  L T.  M BY  PJT  II., 

VI  I • AN  MITT  IMG  Tt  L F 

SELECTION  PROCFL3UPF  IN  SOLCCR 


, H - : T I C C T.Jj  ,iIE  : PS  1 , VT  = >,  M 
It  t ' L i ■ I..  f LLlPfl 


V b K 1 . < 

ZfcRC  I 


FUNCTIONS  <E . 

armoni; 


1 ofl 


? 3?  W J •*  ) - 1 '0  0 
N ’Of  = N3RF  - 1 
2 4 C C J M I I N J - 


se  . r k.  ; 2. 


I IP  DATA  FOR*  AN!  ASSIGN,  I 1CTI 
ru  Ct  uUV.NS  UF  'AT'RIX  AS  I I TAbLL.  1 -VERTlX 

FU  .CT  IONS. 


24  I f ( IDU  >37,37,23 
2 3 NNN  = N3RI  /N3S 
NN  = N 3<F  -NN.,*NSS 
. ' 3 3 J = 1 , . S S 
MAC'  I J ) = NNN 
I r ( J -MM  ) 2'T,  2 i,  3 

2 i NADI  J > = 4\l N -•  1 

3 cdntinj: 

DC  AD  J = 1 , .S3 

NADI J»  * NO.  OF  VERTEX  FUNCTIONS  . J I ■■  FROM  V TEX  J 


< X = N A D ( J ) 

! A D < = 1 , K 
N i-1  = N ’ + 1 
WE  X = N V 2 R I J ) 

A L F A = A N A J ) 

: ; = K * P I / A L : A 

:-r  = t 1/  ( 2*ALFA) 

DO  T 3 { 31  , 32  03 , 34)  , w..  X 
31  CA=W 
CB  = 1 


30  T 3 3 5 
3?  C A = 0 


C E = E F 
CO  TO  35 
3 3 C A = 1 


: e = l = 

30  T 3 3 5 
3 A 3 A = 1 
C 8 = 3 
CE  = EE 

: 21’:  l 1 » M s ) : A 
CDF F ( 2 , J 3 > -CD 
N A H > > = J 
; ( 'P  >=: 

N F U I 'J  3 ) - 4 


:-i  ;k  i dr  : jteg  values  or  L vmi  a=c  , ini  : . icati  * v 

INCREASING  STORAGE  VALUE  E(NP)  LY  K ^ . 

nx=;:+:.do  coci 

I F ( C " - N X - C . i 0 31  >36,36,4 
36  _('>»>=::+]  odd 
40  CON  T 1 N J _ 

s„t  up  columm  lp  fur  r.h.s.  or  ecu at i on . 

_ 7 = J 5 ♦ 1 
N F 0 I L P ) = 3 
N A ( . > ) = 1 : 

3 7 N:1=C 


* . ’ * •* 


SECTION  3.  ' . ! JP  COLUMN  ETS  HP,  EACH  OMPRISING  Til 

LS.  , MATCHING  TwO  rys.  F OR  HARM. WIG  SETS  Ml 

O' J S 1 OF  J . 


SAN 

DO  i>  0 M ! 

650 

' ) 

651 

MM  = NHSt 

652 

DO  50  3: 

65  3 

N ? = V’  + 1 

6 56 

VI  = MS- 

655 

NA  ( \P  ) = 

6d  6 

IF ( ,1YJ 

65  r 

61 

m r u ( J =»  > 

638 

. ( l3  ! =2 

65  ) 

SO  TO  6 

663 

62 

NFJt  NP) 

66  1 

x,<^  A ! ( 

662 

ir ( xx  )6 

66  3 

6 5 

\ f J ( N 3 ) 

66*t 

66 

E ( N 3 1 = M 

663 

1 F ( N T Y 5 

66  6 

6 5 

i r ( j - 1 ) 

667 

6 3 

: ( VI P 1 = - 

668 

66 

+ 2 

669 

50 

CO NT  I NO 

> } J 

670 

nIITt  (3 

671 

hR  11  l 

672 

A < I T ! (6 

673 

LL-  H*_ 

676 

,]  906 

67  5 

LL-  ' ’ *■  1 

6 76 

< < = 2 Z 

677 

906 

. I 

678 

m R I T E I 5 

679 

. R I TE  ( £ 

630 

■ : i ( £ 

68  1 

DO  908 

632 

903 

h R I T E « 6 

6F  i 

a .<  I T l ( 6 

63  6 

A X I IF ( 5 

686 

M I TE  I 6 

636 

90  2 

FORMAT  ( 

1 .olmat 

687 

79  9 

FORMAT  ( 

638 

90  3 

FORMAT ( 
ItfriF  1 X 

63  9 

905 

FORMAT ( 

6 9. 

9"  6 

FORMAT ( 

69  1 

90  7 

F IRMATI 

6 92 

90  9 

FORMAT! 

69  3 

91  C 

FORMAT! 
* 6 X 
2 t6Xi 

696 

RlT JRN 

6 96 

E ND 

Full  0 i ; F u,|5' 


>R1NT  JT  01  DcTaILI  - NFP  FIX  P Pi  I ;TS  AMO  LP  CMlJ":. 


0 IN 


' Xf-  I X 
3 i 3 . * ) 


I ' , / 


MNP) 


Jl  L-UUATI  MS  IN  HAMM'  NIC  S 1 T JH  =*,!,/ 


- f ~ - -”33-: — 


■4&*' 


C 

L L'v 

C 

VA  AM 

c 

6 96 

SJ3  oar  I 11  [ 

6)7 

X X < = ( X<-XJ_ 

69b 

Y Y<  = - ( XK-  X J. 

699 

SSS  = 3X<-3  J_ 

700 

L-'"  = -X7  (-E  ’ 

731 

I F ( M r - ? ) 1 , , 

702 

4 

9a=:yy*ginc 

733 

vr  = f'lM*GJSt. 

734 

jJ  T 3 3 

705 

1 

i 

V A = E MM* 3 2 1 

7C6 

V2  - I . 

72  7 

SO  1J  3 

738 

O 

c 

VA=-EM*EMV  ' : 

739 

V D = ? v|  * w ' v 

7 13 

3 

R-T JIM 

711 

L 'ID 

c 

r 

o 

c 

712 

SUBR3UT I Ml 

713 

XA=  X<  - X J E 

714 

Y A = Y < - Y J E 

715 

Call  gp  jl  a 

r 

-i 

716 

8 A = 8 A - 3 J _ 

71  7 

DDEL= 13*3EL 

7 lb 

I F ( 94+0  0 lL) 

c 

719 

23 

F A = t.  A ♦ 2 * P I 

~T  p 

2 ) 

GS=3XK-3JE 

721 

I F ( 5 Y-202  3 ) 

c 

722 

3 J 

9 A=  34 

723 

1=1 MT-2)5  . . 

7?  4 

4 1 

9A=-AL3G(R- 

725 

GC  TO  50 

726 

3 1 

9 A=C3S( SS-u 

72  7 

GO  TO  53 

728 

£ 

32 

M 5 S = 1 + : M/  1 

r 

c 

729 

EMsSM-I3 

733 

xx=: 

731 

DO  10  <=  1 , . 

732 

1^=1^ 

733 

I F ( M > S-  ’ ) 4 

7 34 

4 4 

t M=  lf|  ( 

735 

45 

, my  = : r»-  ir  ♦ ; 

F.  FUMCTI  MS  t EC 

fi  ST  D M OtFFS.  Of  JMKIMGV  FU  )CTI<  . 

, F( XKtYK, BXK» XJE ,YJEtI  JE,EM, PI,MTf  VA, V > 
fZ,  -I  oJL  ) 4 ( YK- YJT  ) *SI  .(  J!  ) 


>OLW.  i/ERl  < function; 

CCEFF.  (50)  01  . A ' . 


T A3LL 


1,9. 


F i YJ  i’lJL AMGLC  JA>  — ii,  CL 


f JMC  T I . 


r T A I NDI  CAT  Cl'  ' ’ >220: 


SOT  BY  ICAL  DIFF.  INTE RV  '* L D(  L l A 

GIVES  VERTEX  FUNCTI  V S - (HI  |f»!  - I . 


« 


> v • . . A i to. 


mm 


5 0 

ODD' 

. r 

. ' ol 

.6  5:  ' 

- c ^ 

,(.6L  r - 4 

51 

0. 

00  001. 

.0  DOE  .C 

. 

oc 

. 1<  > 

52 

• 

0 23  i 

' • 0 C 1'  0 L - 6 

• 

06 

. 1 

63 

• 

. 0 0 0 C C 

00 

o.oo  ol.  ■: 

. 

.4  :jlt  10-01, 

MT  = 1,2.0R  3 J 


• J A S I 0 R r S 


[ , L..  LT  a ) 


1F(VT-3>33»36»3G 
EMV==M-MT  ♦ - 

c = 0 

I F ( ? < - J . 0 0 '01)20. 20, 1> 

< E = < <** .MM 

1F(MT-2)1,  ,3 

xx  = xx+i  i-ik*k  ) *(  CA*sr.  i • ;-i  * a )-cb*cos( f m*»  a)  ) *»  c 
30  1 3 ID 

XX=XX  + ((-l)**K)*(CA*C:ij(rM*t  3 ) + C3*5  I \ (EM*n 'A  ) ) *7  c 
j3  TD  1 : 

XX  = X X f ( { - 1 ) * * K ) * ( -C A * S I . < ( ! v ‘ *0/  + S S ) ♦ CD  *C  OS  ( f MM#BA+S5  ) ) * E . * u i 

: cm i ^ j: 

\/  a = X X 

< E 1 J \ N 
t xi  P 


P.-.AP  MAPPI  10  Z =C  COSH!,  ) CURVED  ED3E-FJ  'JCT  IONS  SEC. 

V A M i)  V STL)  EF  F S.  CF  UN  t . I N EQ  ( 89 J F OR  MI  = 1 * 2 , 


UBnjriNE  PIYAP  (XK,  YK  ,D  XK  , XJE  ,Y  JE  , J , E M , **  T , JE  , A , \!  1) 

4 X( 11) ,Y( 111 ,C( 10,6) ,A( 1C) ,B(2,1C) ,ANG(LO) ,E< 13  ),  MS(A,1  ) 

C3(  > 0,  1 0 0)  t ML  ( A , 100  ) t At  ( 1 1 ) ( 1 1 ) , X (U),YE(11),XFIX(5.  ) , Yf IX ( 53 

L ( , iO  ),  YEL  ( 2,  50)  » EL  ( 2 , 50  ) , PAft<  2 , II  ) Dl  , )FI  ;t  ) , T . ( a ) 


C(  51  ) , NBDY » NFPOL 


75  3 

COMMON  3F I ( 50)  ,Z fc  < i 1 ) , 

75  A 

::j  'MD  i (5  ) , ICLD(  A ) , . 

f bt> 

: OM  D , H 1111,5:5  1 

756 

; MV  dm  )DK  ( 1 1 ) , 'B  ( 1 1 ) , 

757 

: JT'J'J  MFP  ,\l  ^,5  NSS,  .6 

758 

XX=X<-XJ5 

750 

Y Y=  Y< - Y JE 

763 

2 < = 0 J E 

751 

v[)  3 = 0 

762 

VAA  = 3 

C 

E * U A T 1*1  \i  j ( 

753 

XA  = XX^E)Sl  <)+YY*SI  J(ZK 

76  A 

YA=-XX*SIN(ZK)*YY*CUS(Z 

765 

C:  = AE ( J > * 2-BE ( J1  )**; 

766 

C 1=CC-  <A*<  ?-YA*>*2  ) 

767 

r=(-:i  + sqp 1 1 ci * *2  +a*cc  * 

76  f 

IF!  A 3 S ( Y A )-G.l_-06)A6C, 

76  0 

A3  3 

i = -i«-xa*xa/cc 

773 

IF  < 1)  A 3 1 » A.  1 , A8  2 

77  1 

A 3 1 

r = o 

772 

A 8 7 

T 1 = S D3T  ( r )+S(JR  Ill+ri 

773 

X I * A _ 3 s ( T 1 ) 

77  A 

XA=  XA  *(  3 X-M  2*X1  )-  1 ) 

775 

Y . = Y A * ( EXP(2*XI  ) + 1 ) 

776 

CALL  C'PJLA  ( XA»  YA,RR,  1 A 

777 

E T A * B A 

C 

st  t up  E^smiu 

7 76 

V A * 0 

770 

;p.  = r 

730 

MCY=1 

731 

I F ( J E - 1 ) 1 , . , 3 

7P2 

1 

X 

II  ( ;TY*Ml)-3)3»2t3 

783 

2 

NCY  = ? 

78  A 

3 

03  MO  N*1.NCY 

( y \**2 ) ) ) / ( 2*CC ) 
A r 3 , A 8 2 


* • ’ -V  * •’ 


*S£* 


i. 


785 

IN  '3-  9 99  ) < , 6,  6 90 

736 

690 

ir  ( v, r - ? )6  91 , 502  .6  92 

73  7 

6 9 1 

VAA  = .\L  33  ( <1  ) 

738 

3D  TD  6 06 

73  1 

6 9? 

V 44  = 3A 

7)0 

DU  TD  » 06 

791 

6 

If  ( » - 2 ) i » 6 , 6 

792 

6 

79  3 

!) 

L T A V = p A v o M 

7)6 

X I X = X I * 6 M 

796 

1 \=E  XP(  <1 ' ) 

7 96 

F3  = 0DS(ETA-  ) 

797 

r:=SI\i(fTA  ) 

798 

IF ( MT-2 ) 501,502*503 

7 99 

5 0 1 

V 4A=FA*F3 

8.  DC 

YFa=-F4*rC 

o r 1 

jO;  TD  506 

802 

502 

X l = s I\l( FT  A ) 

6 0 3 

X 2=  C DO  ( L T A ) 

836 

X 3=  XP ( XI  ) 

806 

X 6 = 1 / X 3 

8C6 

Xb  = X3-X  4 

60  7 

X6= X 3 + X 4 

338 

DX=  0. 5*  ( X6< X6-6* X2 *X2 

839 

I F { 43L(  jy.  ) — 0 • 1C  -07)  50 

810 

503 

rx=  .10-07 

811 

50b 

0u:  = 53RT  (OC) 

812 

Ul=X5*X2/ ( CD*L)X ) 

8 1 3 

; 2=-X<  0=X1/(CO*DX) 

8 16 

I F ( l 3-  999  ) 6 94,4  34,4  9 _j 

816 

6 ) _s 

X2=0)S[  3XK  ) 

316 

X 1 = s n ( X K ) 

6 17 

X 3=  X 1 * X I +1-  ! A*C  T , 

818 

X 4=  X I / X 3 

819 

X 6= : TA/  X3 

820 

V A A = - X 1 * ( X4*D1+X5*D? ) 

6 21 

GO  ID  5 D 4 

8 22 

6 96 

FQJ=?X<-8J . 

823 

F 8 = 0 D S ( ETA  -FQJ ) 

82  6 

f:  = sin( eta  -fuj ) 

823 

r A=  F A*E 8 

826 

VAA=FA<t(-Pl*FC+l'2<‘Fi-!) 

8 2 7 

YB?=*-A*  ( -Dl*Fu- 02  *FC  > 

323 

^ TD  304 

329 

50  3 

YAA  = " A*  FO 

830 

VB3=:A*FB 

8 31 

6 0 6 

V A=  V A ♦ Y A A 

832 

6 10 

V.  =V3+Y33 

e 3 3 

KET JPN 

e 36 

t ID 

5 , 5 0 6 


:(  - K4*L>2  + X5*i  1 ) 


L 

L 

C 


GI V , POL R R . . i 

I ' J ( * 2 * ° I ) • 


, 8 A ) \ 


OUT  n 3 POL  • (XA,YA,  »BA»PI) 
N=S3R  r ( X a : *2  + YA**2 ) 

ir  ( a 3 s ( <a  )-o.oc  :od  r,  . 

7 XA=XAO.OC  01 


*.<  *11  Vu  Ju  S 


3 3 7 

Q 

o 

i A - A T A \|  ( Y A / X A ) 

BAG 

I F ( X A ) 9 , 9 , K 

S A 1 

7 

3 A * 3 A ♦ P I 

FA2 

GJ  T 3 12 

1 A3 

1 0 

1 T ( Y A ) 1 1 , 1 2 , I 2 

BAA 

1 1 

*1  A = :•  A ♦ 2 # 5 I 

8 A 5 

1 2 

r fc  r j r m 

8A6 
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A V‘  0 (1)  P A R A f ' :1  r RS  °AR(1»J)A.D  PAR(?,J)  A 40 

SLOPE  S ! ( 1 , J ) A JD  8 (2 . J)  AT  E ND  POINTS  f 1 Y. 

3 ML  vi  I r OR  KK=  ; 

(2)  30i?R3S.  (XK,YH  AND  S L PL  BX  K AT  KK 
' INTS  ON  SEGMENT, KK>0.  STORE!  I . <EL  » Y L»BEL. 


3 A 7 

SJR  OUT  IME 
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S A 3 

0 JMM 3R 

X(  1 1 

) , Y( 

CG(  79,130)  , 

J L ( A 

C , X E L ( 2 , 
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•YlLI 

3 A 7 

13 r I ; 

( 50  ) 

S50 

: 7"  r. 

( 5 

) , A rl 

° 5 1 

3 3‘) v 3'4 

4 R G . 

111) 

3 5? 

OO.MMDM 

M 3 K ( 

11), 

8 5 3 

3 3 M M 3 '4 

4”  n , 

•■4  L Z R 

8 5 A 

1 F ( K < ) 3 

03  , 

■33  , 3 

L ( 2 ,5 0 ) » P AR  ( 2 » 1 C ) *GMOD*  Of  1 < ( 


( I 1 ) , JVE  < ( 1 1 ) , ..  (A  ) , IP,  •;  AX  I S,  vi 
\F  S( 1 1 ) , O Y ( 1 1 ) , TYP(  1 1 ) , JF  IX  <51 
, 13  tL  ( 1 . ) ,L3C  (53  ) OY.AJFP'.'L 


o ) 

) , YF IX  < 50 
OR  ( A ) 

> , I ( * 0 ) 


‘ I X ( 50  ) 


P AR  AM  L T 

I F ( RTY-M  J 1-3)302,301  ,301 
> A R ( 1 , J ) = C 
•>  A R ( 2 , J ) = 2 5 I 
A ( J ) =?*r,I 

GO  TD  5 DO 

I F ( • 5 ( J > ) 3 3,  3 0 3, 3 G A 
6 E ( J I = A L ( J ) 

00  3)5  K= l , 2 

J J = J + < - 1 

XX  = X (JJ ) - X tJ) 

Y Y = Y ( JJ  )- Y,  ( J ) 

2 < = 7 L ( ) ) 

xA=xx*:os(cX)  + YY*sri(ZK) 

Y A = - X X * S I M Z < ) + Y Y #C0  S 1 Z P ) 

X XA  = 3 l ( J ) *XA 
YYA=  AL(  J ) *YA 

;all  :5  ol  a ( xxa  ,yya , r.’ , l a,  pi  » 
XA=  C3S( 5 A ) 

Y A = S I M ( 3 A ) 

IF  ( YA  ) A 33  , a 32 ,4  .3 
YA  = YAO.OO  JLl 
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